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When solving numerical constraints such as nonlinear equations and inequalities, solvers often 
exploit pruning techniques, which remove redundant value combinations from the domains of 
variables, at pruning steps. To find the complete solution set, most of these solvers alternate the 
pruning steps with branching steps, which split each problem into subproblems. This forms the 
so-called branch-and-prune framework, well known among the approaches for solving numerical 
constraints. The basic branch-and-prune search strategy that uses domain bisections in place of 
the branching steps is called the bisection search. In general, the bisection search works well in 
case (i) the solutions are isolated, but it can be improved further in case (ii) there are continuums 
of solutions (this often occurs when inequalities are involved). In this paper, we propose a new 
branch-and-prune search strategy along with several variants, which not only allow yielding better 
branching decisions in the latter case, but also work as well as the bisection search does in the 
former case. These new search algorithms enable us to employ various pruning techniques in 
the construction of inner and outer approximations of the solution set. Our experiments show 
that these algorithms speed up the solving process often by one order of magnitude or more 
when solving problems with continuums of solutions, while keeping the same performance as the 
bisection search when the solutions are isolated. 
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1. INTRODUCTION 

A constraint satisfaction problem (CSP) consists of a finite set of constraints spec¬ 
ifying which value combinations from given domains of its variables are admitted. 
It is called a numerical constraint satisfaction problem (NCSP) if the domains are 
continuous. NCSPs such as systems of equations and inequalities arise in many 
industrial applications. A method for solving NCSPs is called complete if every so¬ 
lution can be found by it in the infinite time and can be approximated by it within 
an arbitrarily small positive tolerance after a finite time, provided that the under¬ 
lying arithmetic is exact. Most available complete solution methods are instances 
of the branch-and-prune framework, which interleaves branching steps with prun¬ 
ing steps. Roughly speaking, a branching step divides a problem into subproblems 
whose union is equivalent to the initial problem in term of the solution set, and 
a pruning step reduces a problem in some measure. The reader can find a more 
detailed discussion on branch-and-prune methods in [Vu 2005, Section 3.2]. 

The need for completeness arises in many applications such as safety verification 
and computer-assisted proofs [Schichl 2003, Section 3.1]. In design applications such 
as estimation and robust control [Asarin et al. 2000; Jaulin et al. 2001], automation 
and robotics [Jaulin et al. 2001; Lee et al. 2002; Lee and Mavroidis 2002, 2004; 
Neumaier and Merlet 2002], civil engineering [Lottaz 2000; Vu 2005], and shape 
design [Snyder 1992], the solution set of an NCSP often expresses equally relevant 
choices that need to be identified as precisely and as completely as possible. In 
a design process, one often desires to find as many solutions as possible because 
investigating many solutions at earlier stages potentially increases the chance of 
success at later stages. This also allows identifying good design choices. Thus, 
the complete solution set is often sought for, provided that the response time is 
reasonable. Since a general NCSP is NP-hard, the time for finding all solutions 
at a high precision is often prohibitively long. In most cases, a low or medium 
precision is however sufficient for applications. Hence, there is a tradeoff between 
timely but less precise information and slow but more precise information. 

The necessary background is presented in Section 2. The rest of the paper is or¬ 
ganized as follows. We first prepare general settings about solution representations 
and formal definitions of domain reduction and splitting operators in Section 2.4 
and Section 3, respectively. We then present, in Section 4, a generic branch-and- 
prune search algorithm, called BnPSearch, that enables the incorporation of do¬ 
main reduction and splitting operators (Section 4.1), and then present, as instances 
of BnPSearch, the bisection search (Section 4.2) and two new search algorithms, 
called UCA5 and UCA6 (Section 4.3). Roughly speaking, the effectiveness of the 
new algorithms is mainly due to the following policies, if some constraint, C, is 
predicted to contain continuums of feasible points: 

—Working on the negation of constraints to find feasible regions. Apply 
domain reduction techniques to the negation of C. Let x be the input domains 
and x' the resulting domains. Then x* := x \ x' is feasible w.r.t. C. Thus, C 
can be immediately removed from subproblems defined on subregions of x*. 

—Reducing the number of variables in subproblems. In a subproblem gen¬ 
erated during the solving process, some variables may not occur in any constraint 
under consideration and thus no longer need to be considered. 
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Taking the influence of all constraints on each other into account at 

each iteration. Potentially, this reduces the search space better than solving 

each constraint at a time does. 

A simple heuristic to predict whether a constraint contains continuums of solu¬ 
tions is to check if this constraint is an inequality. In general, UCA5 and UCA6 
allow better branching decisions than the basic domain bisection. Although pro¬ 
viding an accurate representation of solutions, in some cases these algorithms are 
still slow and provide verbose representations. The first reason is that the orthog¬ 
onal splitting policy in these algorithms generates a significant number of nearly 
aligned boxes near the boundaries of constraints. The second reason is that these 
algorithms often have to spend too much effort on producing too small boxes with 
respect to the precision e, which is predetermined by users. 

We later propose, in Section 5, an improved instance, called UCA6+, of 
BnPSearch to tackle the above two limitations. The improvement in UCA6^ is 
twofold. First, UCA6^ utilizes domain reduction techniques better when the preci¬ 
sion is recognized as being sufficient. Namely, it tells domain reduction techniques 
to avoid unnecessarily spending too much effort on reducing the domains whose 
sizes are smaller than e. When the sizes of a certain number of domains are smaller 
than e, UCA6+ allows resorting to a simple solver that is more efficient for small 
and very low dimensional problems. The gain is then in both the computational 
time and the alignment of boxes. Second, UCA6+ allows resorting to geometric 
representation techniques to combine aligned boxes produced in the previous stage 
into larger equivalent boxes. The representation of the solution set is therefore 
more concise. This potentially accelerates querying and complicated operations on 
the explicit representation of solutions. 

In general, our new search algorithms improve the solving process when there 
are continuums of solutions, and keep the same procedure and performance as the 
bisection search when the solutions are isolated. In the former case, our search al¬ 
gorithms allow producing inner and outer approximations w.r.t. a predetermined 
precision e. Moreover, a large percentage of the outcome are often proved to be 
sound solutions. Our experiments in Section 6 show that the new search algorithms 
significantly improve the efficiency as well as the conciseness of solution represen¬ 
tations. The conclusion is finally given in Section 7. 

2. BACKGROUND AND DEFINITION 
2.1 Numerical Constraint Satisfaction 

We recall in this section two central concepts of constraint programming. 

Definition 2.1. A constraint on a finite sequence of variables, (a;i,... ,Xk), tak¬ 
ing their values in respective domains, (Hi,... ,Dk), is a subset of the Cartesian 
product Di X ■■■ X D^, where /c € N (N is the set of natural numbers). 

Definition 2.2 CSP. A constraint satisfaction problem, abbreviated to CSP, is a 
triple {V,D,C) in which V is a finite sequence of variables (fi,..., u„), 22 is a finite 
sequence of the respective domains of the variables (ui,... ,u„), and C is a finite 
set of constraints, each on a subsequence of V. A solution of this problem is an 
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assignment of values from T) to V respectively such that all constraints in C are 
satisfied. The set of all solutions is called the solution set. 

The reader can find many more concepts in [Apt 2003]. In this paper, we only 
focus on numerical CSPs, which are defined as follows. 

Definition 2.3 NCSP. A domain is said to be continuous if it is a real interval. 
A numerical constraint is a constraint on a sequence of variables whose domains 
are continuous. If all constraints of a CSP are numerical, it is called a numerical 
constraint satisfaction problem (abbreviated to NCSP). 

An NCSP can be viewed as a constrained optimization problem with a constant 
objective function. Thus, it can be theoretically solved by using mathematical pro¬ 
gramming (MP) methods for solving constrained optimization problems. However, 
most of the efficient MP methods are heavily based on the influence between ob¬ 
jective functions and constraints, thus not efficient for directly solving NCSPs. 

Since thirty years ago, constraint satisfaction techniques have been being devised 
to solve CSPs with discrete domains. These techniques perform reasoning proce¬ 
dures on constraints and explore the search space by intelligently enumerating so¬ 
lutions. In order to solve NCSPs by means of constraint satisfaction, continuous 
domains have often been converted into discrete domains by using progressive dis¬ 
cretization techniques [Sam-Haroud 1995; Lottaz 2000]. Still, these methods are 
often inefficient. Later on, many mathematical computing techniques for contin¬ 
uous domains have been integrated into the framework of constraint satisfaction 
in order to solve NCSPs more efficiently. Nowadays, these techniques are often 
referred to as constraint programming techniques, which imply the combination of 
computing and reasoning aspects. A much more extensive discussion on numerical 
constraint solving can be found in [Vu 2005]. 


2.2 Interval Arithmetic 

Let K-oo = K. U {—oo,-|-oo}. The lower bound of a real interval x is defined as 
inf(x) e K-oo) and the upper bound of x is defined as sup(x) € Moo- Let denote 
X = inf(x) and x = sup(x). There are four possible intervals x with these bounds: 


• The closed interval defined as x 

• The open interval defined as x 

• The left-open interval defined as x 

• The right-open interval defined as x 


[x,x] 

]x,x[ 

]x,x] 

[x,x[ 


{x G M I x < a; < x}; 
{x G M I x < X < x}; 
{x G M I X < X < x}; 
{x G M I X < X < x}. 


Let I be the set of all closed intervals and lo the set of all intervals. The interval 
hull of a subset S' of M is the smallest interval (w.r.t. the set inclusion), denoted 
as ns, that contains S. For example, n(]l,3] U {2,4}) = ]1,4]. Given a nonempty 
interval x, we define that the midpoint of x is mid(x) = (inf(x) -|-sup(x))/2 and the 
width of X is w(x) = sup(x) — inf(x). We also agree that w(0) = 0 and mid(0) = 0. 
A box is the Cartesian product of a finite number of intervals. The concepts of the 
midpoint and width are defined on boxes in a componentwise manner. 

Fundamentally, if x and y are two (real) intervals, then the four elementary 
operations for idealized interval arithmetic obey the rule 


xoy = {xoy I X G X,?/ G y}, Vo G {-k,G-}. 
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Thus, the results of the four elementary operations in interval arithmetic are exactly 
the ranges of their real-valued counterparts. Although the rule (1) characterizes 
these operations mathematically, the usefulness of interval arithmetic is due to the 


operational definitions based on interval bounds. For example, let x = [x, a;] and 
y = [t/,y] be two closed intervals, interval arithmetic shows that: 

x-Py = [x + y,x + y]-, ( 2 ) 

x-y = [x-y,x-y\] ( 3 ) 

X * y = [min{xy, 2 y,Ty,^},max{x|;,xy,Ty,xy}]; (4) 

X -P y = X * (1/y) if 0 ^ y, where 1/y = [l/y, l/y\. (5) 


Simple arithmetic expressions are composed of these four elementary operations. 

An interval form, f : I™ ^ I", of a real function f : D C ^ R" is constructed 
conforming to the inclusion property, the value of the interval form encloses the 
exact range of the real function, that is, Vx C I™ : f(x) C f(x) or, equivalently, 
Vx G I™, X G : a; G X f{x) G f(x). 

The finite nature of computers precludes an exact representation of the real num¬ 
bers. The real set R is therefore approximated by a finite set F of floating-point 
numbers [Goldberg 1991], including —oo and -Poo. The set of real intervals is then 
replaced with the set !<> of closed floating-point intervals with bounds in F. The 
interval concepts are similarly defined on while conforming to the inclusion prop¬ 
erty. The power of interval arithmetic lies in its implementation on computers. In 
particular, outwardly rounded interval arithmetic allows computing rigorous enclo¬ 
sures of the ranges of functions. An interval is said to be canonical iff it does not 
contain two different intervals whose union is not an interval. A box is said to be 
canonical iff all its intervals are canonical. 

The reader can find extended introductions to interval analysis in [Moore 1966, 
1979; Alefeld and Herzberger 1983], interval methods for systems of equations in 
[Neumaier 1990], interval methods for optimization in ]Hansen and Walster 2004], 
and some recent applications of interval arithmetic in [Jaulin et al. 2001]. 

Most interval arithmetic libraries have been implemented for closed intervals only. 
One however can use these libraries to perform some computations on open/closed 
intervals. Indeed, every interval x with bound values x and x is contained in the 
corresponding closed interval [x, x] G I. If the computations are domain reduction 
or complementary boxing, we can perform these computations on the correspond¬ 
ing closed intervals of the domains to get new domains, and then take the set 
intersection of these new domains with the initial general intervals. For example, 
after performing a domain reduction technique on the closed interval [s, T], we get 
a closed interval [y,y] C [x,x]. The result to be obtained is the set intersection 
X n [y,y], thus may be open or closed. 

2.3 A Short Overview of Branch-and-Prune Solution Methods 

To be able to find the complete solution set of an NCSP, most solvers follow the 
branch-and-prune framework, a well known approach for solving numerical con¬ 
straints. In this framework, a solver alternates pruning steps with branching steps 
until reaching the required precision, where a pruning step attempts to reduce each 
considered problem in some measure and a branching step splits each considered 
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problem into subproblems. The basic branch-and-prune search that uses domain 
bisections in place of the branching steps is called the bisection search. 

A pruning technique that attempts to reduce the domains of problems without 
discarding any solution is called a domain reduction technique. In contrast to a 
domain reduction technique, a test (e.g., an existence test, a uniqueness test, an 
exclusion test and an inclusion test) does not change the domains of an input 
problem but maps this problem to a predetermined status. In particular, existence, 
uniqueness and exclusion tests are to check if a problem has at least one solution, 
a unique solution and no solution, respectively. The outcome of these tests is thus 
a Boolean value. Quite differently, an inclusion test is to check if all points under 
consideration are solutions, non-solutions, or else; which is defined as follows. 

Definition 2.4. Let X be a sequence of n real variables. An inclusion test is a 
function t that takes as input a (domain) box x G I" and a finite set C of constraints 
on a subsequence Y oi X (assuming Y is well defined on x) and that returns either 
feasible, infeasible, or unknown such that: 

(1) If r(x,C) = feasible, then every point in x satisfies all constraints in C. 

(2) If r(x,C) = infeasible, then no point in x satisfies all constraints in C. 

The inclusion test t is said to be trivial if it always returns unknown. It is said to 
be e-strong if, for every x, the truth of w(x) < s A r(x,C) = unknown implies the 
existence of a feasible point of C in x. 

An extended overview of fundamental and recent complete methods for solving 
NCSPs has been presented in [Vu 2005, Chapter 3]. In summary, those methods can 
be viewed as instances of the branch-and-prune framework. Many of them integrate 
existence, uniqueness or exclusion tests at pruning steps of the bisection search. At 
a branching step, they bisect a domain, which is often chosen as the largest with 
respect to some measure (e.g., domain size), provided that this domain is amenable 
to be split (e.g, its size is greater than a predetermined precision s). Because 
most solution methods have been designed to solve a square system of equations, 
they only aim at generating a collection of tiny boxes, each encloses a solution. 
This approach is referred to as the point-wise approach. It may be reasonable for 
solving NCSPs with isolated solutions (see Figure la), but are often inefficient, 
when applied in a straightforward manner, for solving NCSPs with continuums of 
solutions (see Figure lb). In the latter case, neither the computational time nor 
the compactness of the solution representation are satisfactory. 

It is possible to enhance the solving process in the point-wise approach by re¬ 
placing the existence, uniqueness and exclusion tests with domain reduction tech¬ 
niques in a straightforward manner. In the rest, the resulting search strategy will 
be called dichotomous maintaining bounds by consistency (DMBC). When solving 
NCSPs with non-isolated solutions, DMBC search techniques may be able to cover 
a spectrum of non-isolated solutions with a number of subsets of M”. However, it 
is often not possible to prove that a subset of the outcome are all solutions. 

In contrast to the point-wise approach, another one, called the set-covering ap¬ 
proach, has been developed in order to represent continuums of solutions more 
reasonably. It aims at covering, as accurately as possible, continuums of solutions 
with inner and outer approximations, each consists of a number of subsets of M". 
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Fig. 1. (a) An NCSP with four isolated solutions (grey dots); (b) An NCSP with two continuums 
of solutions (grey regions). 


All points represented by the inner approximation are proved to be solutions. Usu¬ 
ally, the representation of inner and outer approximations is made simple such that 
the costs of usual operations on this representation are as cheap as possible. The 
subsets in these approximations are often chosen to be simple ones (e.g., boxes). 
Owning to their nice properties, boxes have been used in many set-covering tech¬ 
niques, thus forming the so-called box-covering techniques. 

For simplicity, in this paper we restrict our attention to solution meth¬ 
ods that use boxes as elements of the outcome. One has often employed 
inclusion tests at pruning steps of typical box-covering branch-and-prune methods 
to prove that all points in some subsets are solutions. In this paper, a DMBC-like 
search technique that combines both domain reduction techniques and inclusion 
tests at its pruning steps is called a DMBC+ search technique (see Section 4.2). 
Therefore, DMBC^ search techniques may, depending on the strong of employed 
inclusion tests, be able to provide inner approximations of the solution set. We 
hence say that DMBC+ search techniques belong to the box-covering approach 
and that DMBC search techniques belong to the point-wise approach. Note that 
both DMBC and DMBC^" search strategies do not remove constraints from con¬ 
sideration during the solving process. That is, they always consider the whole set 
of initial constraints when resorting to domain reduction techniques and tests. 

When solving an NCSP with continuums of solutions at a high precision, DMBC 
techniques often provide a huge collection of tiny boxes as an outer approximation of 
the solution set while DMBC+ techniques are usually able to provide both inner and 
outer approximations of the solution set, each approximation is a much more concise 
collection of boxes. However, the number of boxes used by DMBC+ techniques to 
approximate the boundary of continuums of solutions is still very high in many 
cases. Therefore, either their applicability is restricted or the tractability limit is 
rapidly reached. The new search techniques proposed in this paper will tackle this 
limitation (see Section 4 and Section 5). 

Because most computers were equiped with floating-point number systems, an 
important issue is how to deal with rounding errors occurring in computations on 
floating-point numbers. Owning to the inclusion property of outwardly rounded 
interval arithmetic, one has often been using it to tackle the issue of rounding 
errors. Recent interval arithmetic based methods are able to solve a number of 
NCSPs efficiently while still enjoying the completeness. Those methods can be 
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viewed as instances of the branch-and-prune framework, although they mainly focus 
on improving the pruning steps. 

Interval constraint solvers such as CLP(BNR) [Benhamou and Older 1992, 1997], 
Numerica [Van Hentenryck 1998] and ILOG Solver [ILOG 2003] have shown their 
ability to efficiently find all solutions of certain instances of NCSPs within an 
arbitrary positive tolerance. They are instances of the branch-and-prune frame¬ 
work, and most of them use a simple branching policy like the domain bisection 
or dichotomization as their default branching policy while leaving more advanced 
branching policies to users. In other words, they essentially follow the point-wise 
approach and only aim at solving NCSPs with isolated solutions. 

Recently, a number of box-covering methods have been developed in [Jaulin 
1994; Jaulin et al. 2001], [Sam-Haroud and Faltings 1996], [Garloff and Graf 1999], 
[Collavizza et al. 1999] and [Benhamou and Goualard 2000; Benhamou et al. 2004] 
in order to represent continuums of solutions. Although those methods are more 
suitable for dealing with continuums of solutions than the point-wise methods are, 
their branching policies still have at least one of the following limitations: 

They are not complete methods in general. For example, Collavizza et al. 
[1999] have proposed a technique to extend a known feasible box of an inequality 
of the form f{x) < 0 by performing box consistency (a kind of domain reduction) 
on its associated equation, f(x) = 0. Unfortunately, their results (Proposition 1 
and 2 in that paper) do not hold for general constraints. 

They are only designed for special constraints. For example, the technique 
proposed in [Garloff and Graf 1999] uses Bernstein polynomials to construct 
algebraic inclusion tests for use in a DMBC/DMBC+ search, and is restricted 
to polynomial constraints. The technique proposed in [Sam-Haroud and Faltings 
1996; Lottaz 2000] is restricted to the class of NCSPs with convexity properties. 
The technique proposed in [Benhamou and Goualard 2000; Benhamou et al. 2004] 
is originally designed for solving universally quantified constraints. 

They do not fully exploit the power of domain reduction techniques. 
Namely, they only interleave inclusion tests with uniformly splitting policies on 
all variables: each box produced by the splitting is tested for inclusion. The 
outcome can be structured into the form of a 2^-tree. In [Jaulin 1994], this 
process is performed in the space of all variables. However, in [Sam-Haroud and 
Faltings 1996; Lottaz 2000], only binary and ternary constraints obtained by 
ternarizing^ the initial NCSP are considered for the construction of 2*-trees. 
They do not sufficiently take the influence of constraints on each other 
into account during the solving process. For example, the solution method 
in [Sam-Haroud and Faltings 1996; Lottaz 2000] construct quadtrees and oetrees 
for individual binary and ternary constraints, respectively, and finally perform 
a constraint propagation on those trees. On the other hand, at each iteration 
the solution method in [Benhamou and Goualard 2000; Benhamou et al. 2004] 
considers a constraint and solves it within every (domain) box produced as the 
outcome of the previous iteration. 

^Ternarizing an NCSP is to recursively replace each binary arithmetic subexpression with an 
auxiliary variable until the arity of all the resulting expressions is at most three. 
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Most methods represent the solution set of an NCSP explicitly in the space of 
initial variables, thus suffering from the high space complexity when there are con¬ 
tinuums of solutions. The only one exception we know of is the work of Sam-Haroud 
and Faltings [1996] and Lottaz [2000], in which the authors have proposed to replace 
the explicit representation (in the space of initial variables) with a semi-explicit 
one, which is maintained by a number of quadtrees and octrees. Although that 
semi-explicit representation reduces the space complexity, it increases the querying 
time for a solution. Notice that the uniformity of splitting must be maintained in 
those methods in order to do propagation among 2^-trees. Therefore, the power of 
domain reduction cannot be fully exploited in this method. 

One of the most recent improvements to search is the work in [Benhamou and 
Goualard 2000; Benhamou et al. 2004], which is summarized as follows. In order 
to find feasible regions of universally quantified constraints of the form \/t £ Dt 
f{x,t) < 0, Benhamou and Goualard [2000]; Benhamou et al. [2004] have proposed 
to perform a kind of domain reduction on their negation, /(x, t) > 0. This operation 
is called a negation test. It encloses all possibly infeasible regions and the remaining 
is feasible. Now consider a subproblem living in a domain box x. A procedure, 
called ICAbSc in [Benhamou and Goualard 2000], takes as input a constraint, C. 
The procedure ICAbSc performs a negation test on C to reduce x to a new domain 
x'. If x' is an empty set, then every point in x satisfies C; otherwise, split x\x' into 
boxes and then dichotomize x'. Now, C can be removed from all subproblems that 
do not have domains in x'. The procedure ICAbSc recursively performs the above 
operations on all resulting subproblems that still have C as a running constraint^ 
and that have domains larger than a predetermined precision e. A search method, 
called ICAb5 in [Benhamou and Goualard 2000], repeats the procedure ICAbSc for 
each constraint, one by one, until all constraints have been processed. See also 
Section 3.2 for more discussion on the negation-based approach. 

2.4 Representation of Non-isolated Solutions 

2.4.1 Inner and Outer Approximations. In case the solution set of an NGSP 
is empty or consists of isolated points, its representation is usually simple. The 
representation of the solution set is not simple in other cases, especially when the 
solution set contains continuums of solutions. In general, the solution set of an 
NGSP is a relation on M", where n is the number of variables in the NCSP. A 
relation can be theoretically approximated by a superset and/or a subset. 

Definition 2.5 Inner Approximation. Given a relation, S C R”, a set S~ C M" 
is called an inner approximation of S if it is contained in S'; that is, S~ C S. 

Definition 2.6 Outer Approximation. Given a relation, S C R”, a set S’*" C R" 
is called an outer approximation of S if it contains S; that is, S’*" A S. 

When a relation on R”, such as the solution set of an NCSP, is approximated 
by an inner approximation and/or an outer approximation. The latter is a sound 
approximation (i.e., it only contains solutions), but may lose some solutions. Con¬ 
versely, the former is a complete approximation (i.e., it contains all solutions), but 


running constraint is a constraint that is currently still under consideration. 
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Fig. 2. An example of inner/outer/boundary union approximations of a circle with interior: the 
collection of the dark grey boxes is an inner union approximation (EB^); the collection of the light 
grey boxes is a boundary union approximation (B3®); the collection of the light and dark grey 
boxes is an outer union approximation (ES®). 

may contain some points that are not solutions. Given an exact representation TZ, 
such as a collection of boxes or a tree of boxes, of a relation S C R", we denote by 
pts(7^) the set of points in S. 

One often uses the volume difference, vol(5'+) \ vol(S'“), to measure the degree 
of mismatch between inner and outer approximations. The exact approximation 
errors are then bounded by this measure. 

2.4.2 Union Approximations. Since the time for querying a point in a box is con¬ 
stant, one often approximates a relation S C M” by a collection of pairwise disjoint 
boxes, where two boxes are said to be (strictly)^ disjoint if they have no common 
points. Such a collection is called a collection of disjoint boxes for short. The rep¬ 
resentation of a collection of disjoint boxes which is constructed by enumerating 
these boxes and storing their coordinates is called the disjoint box representation 
[Aguilera 1998]. Among the approximations by collections of boxes, the following 
three attract the most attention in practice because of their simplicity. 

Definition 2.7. Given a relation S C R". An inner union approximation of S, 
denoted by ffl^[S'], is a collection of disjoint boxes in I" such that S A pts(ffl^[S']). 

Definition 2.8. Given a relation S C R". An outer union approximation of S, 
denoted by ffl‘^[S'], is a collection of disjoint boxes in I" such that S C pts(ffl®[5']). 

Definition 2.9. Given a relation S C R". A boundary union approximation, 
denoted by of S (with respect to an inner union approximation and 

an outer union approximation is a collection of disjoint boxes in I" such 

that pts(ffl®[S']) = pts(ffl®[5']) \ pts(ffl^[S']). 

Remark 2.10. Notice that is not a function, where X G {I,0,Bj. In this 
paper, we will always refer to ffl®[S'] with respect to some ffl^[S'] and some 
even when not mentioned explicitly. 


®The main results in this paper still hold if we relax the disjointness such that disjoint boxes may 
have common points on their facets but not in their interiors. 
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The concepts of union approximations are depicted in Figure 2. Note that we al¬ 
ways have the identity pts(ffl^[S'])npts(ffl®[S']) = 0. In practice, one often computes 
and hrst, and then obtains a^[S] simply by a^[S] := U 

The worst-case query time of a bounding-box tree in is + k), where 

N is the number of boxes and k is the number of boxes intersecting the query 
range [Agarwal et al. 2001]. It is therefore useful to construct inner and/or outer 
union approximations of an unknown relation (e.g., the solution set of an NCSP) 
in the form of a bounding-box tree. That is, the box represented by any node of 
the tree contains the box represented by its child node, and all boxes represented 
by the children of any node are disjoint. Fortunately, this property is enjoyed by 
branch-and-prune solution methods assuming that the domains are intervals. 

Several authors have recently addressed the construction of inner and/or outer 
union approximations of the solution set of an NCSP. In [Jaulin 1994], union ap¬ 
proximations are hierarchically constructed in the form of a 2^-tree in the space 
of initial variables. This technique has shown its practical usefulness in robotics, 
automation and robust control. The method in [Sam-Haroud and Faltings 1996; 
Lottaz 2000] also aims at the construction of 2^-trees in a similar way. However, 
only binary and ternary constraints obtained by ternarizing the initial NCSP are 
considered for the construction. That is, only quadtrees and octrees are constructed. 
The solution set is Hnally approximated by a number of quadtrees and octrees rather 
than a single 2^-tree. The space complexity of approximations is thus reduced. The 
approach is however restricted to the class of NCSPs with convexity properties. 

Most recently, the method proposed in [Benhamou and Goualard 2000; Ben- 
hamou et al. 2004] corrects and extends the idea in [Collavizza et al. 1999] to con¬ 
struct inner union approximations of universally quantified constraints. Namely, 
when solving a universally quantified constraint of the form Wt € Dt : f{x, t) < 0, 
the application of a domain reduction technique on the constraint f{x, t) > 0 allows 
finding feasible regions on which an efficient search is based (see Section 3.2). 

2.4.3 The Precision and Accuracy of Union Approximations. The cost for 
achieving a given accuracy of approximations is often very high. Alternatively, 
most constraint solvers stop splitting a box, which represents the domains of a sub¬ 
problem, as soon as the size of this box is not greater than a given positive precision 
e (and this box is called an e-bounded box). Some other solvers may attempt to 
apply a pruning technique or a test to £-bounded boxes before classifying them as 
undiscernible; thus, the name undiscernible box has come out. 

In general, different constraint solvers use different criteria for leaving £-bounded 
boxes unprocessed. If a technique that is applied to £-bounded boxes before claim¬ 
ing them as undiscernible is used by a solver, then it can be used for the other 
solvers as well. Therefore, the comparison of search techniques should be based on 
the same criteria of classifying £-bounded boxes as undiscernible. We propose to 
use monotonic inclusion tests defined in the following dehnition for this purpose. 
Let x[F] denote the projection of a set x on the subsequence Y of variables. 


Definition 2.11 Monotonicity. Let use the same notations as in Definition 2.4. 
The inclusion test r is said to be monotonic if, for every box x' and every Hnite set 
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C of constraints on Y such that x[y] C x'[y] C HceC' have 

r(x, C) = unknown r(x',C) = r(x^,C U C^) = unknown. (6) 

Once the monotonicity holds for the value unknown of an inclusion test, it also 
holds for other values as shown below. 

Theorem 2.12. Let use the same notations as in Definition 2.11. If t is a 
monotonic inclusion test, then 

—If t{'x.,C) = feasible holds, then t(x',C) = t(x',CUC') = feasible holds 
for every finite set C of constraints on Y and every nonempty box x! such that 
x'[y] C x[F] C flcec' ^ 

—If t{x,C) = infeasible holds, then t{x! ,C) =t(x',CUC') = infeasible holds 
for every finite set C of constraints on Y and every box x' such that x'lFl C 
x[T] C flcec' C holds. 

Proof. This is easily proved by contradiction based on Definition 2.11. □ 

Based on the concept of a monotonic inclusion test, we dehne the precision of 
union approximations (or a solution algorithm computing them) as follows. 

Definition 2.13. Given an NCSP V = {V,D,C), a precision (vector) £ > 0, 
and a monotonic inclusion test r. A solution algorithm that computes inner and 
boundary union approximations is (and thus those approximations are) said to 
be of the precision s w.r.t. the monotonic inclusion test r if the boundary union 
approximation equals (w.r.t. the set union) to a collection lA of disjoint £-bounded 
or canonical boxes in I" such that 

Vx S hi : t(x,C) = unknown. (7) 

If T is trivial, we say for short that the solution algorithm and the computed approx¬ 
imations are of the precision e. If r is £-strong, we say that the solution algorithm 
and the computed approximations are e-accurate. 

It is easy to see that a solution algorithm is complete if it is £-accurate (i.e., r is 
£-strong) for all sufficiently small £ > 0. 

3. REDUCTION AND SPLITTING OPERATORS FOR EXHAUSTIVE SEARCH 

Our improvements to the classic search strategies (i.e., DMBC and DMBC+) will 
presented in Section 4 and Section 5. In order to present those improvements 
uniformly and concisely, we generalize and modify some previously existing concepts 
in the next four subsections. 

3.1 Domain Reduction Operators 

First, we define the concept of a domain reduction operator as follows. 

Definition 3.1 Domain Reduction Operator, DR. Given a sequence X of n real 
variables associated with domains D. A domain reduction operator DR for numerical 
constraints is a function that takes as input a box x G I" contained in D and a 
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finite set C of constraints on X, and that returns a box in I”, denoted by DR(x,C), 
satisfying the following properties: 

(Contractiveness) DR(x, C)Cx, (8) 

(Correctness) DR(x,C) 3 xC n (9) 

cec 



Fig. 3. A domain reduction (DR) operator is applied to a box x and a constraint set {Ci,C 2 }. 


A domain reduction operator has also been referred to as a narrowing operator, 
contracting operator, or contractor in literature. We adopt the terminology domain 
reduction operator^ because it is mnemonic and the terminology domain reduction 
has been widely accepted in many fields, not only in constraint programming. 

Definition 3.2 Monotonicity. Given a sequence X oi n real variables associated 
with domains D. A domain reduction operator fj, is said to be monotonic if, for 
every set C of constraints on A", we have 

Vx,x'gIo, X C H, x'C I? : X C x' ^(x,C) C /i(x',C). (10) 

In constraint programming, domain reduction operators are usually constructed 
by enforcing either box consistency [Benhamou et al. 1994], hull consistency [Ben- 
hamou and Older 1992, 1997], or fcB-consistency [Lhomme 1993]. Although these 
domain reduction operators enjoy the monotonicity, many domain reduction op¬ 
erators do not enjoy the monotonicity but are still very efficient in practice. The 
concept of a domain reduction operator is depicted in Figure 3. Other examples 
of domain reduction operators are depicted in Figure 5. The following property is 
straightforward but interesting for constraint solving. 

Theorem 3.3. Given a set C of constraints on a sequence of n real variables 
associated with domains V. Suppose x G I" is a box contained inV. If there exists 
a domain reduction operatorDR that maps (x,C) to an empty set (i.e., DR(x,C) = 0 
then C is inconsistent in x; that is, 

X C [where = I? \ C). (11) 

CeC 

Prooe. It follows from the correctness of domain reduction operators (Defini¬ 
tion 3.1) that X n Clcec C = 0- Thus, we have x C -nC. □ 
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3.2 Complementary Boxing Operators 

We recall that in [Benhamou and Goualard 2000; Benhamou et al. 2004], a negation 
test takes as input a universally quantified constraint of the form V< S : f{x, t) < 
0 and performs a kind of domain reduction operator on its negation, f{x, t) > 0. It 
encloses all possibly infeasible regions and the remaining is feasible. 

The negation test has been extended to solve classic numerical constraints by 
Silaghi et al. [2001], in which a negation test is applied to inequalities of the form 
f{x) < 0. The search algorithm in [Silaghi et al. 2001], called UCA6, employs the 
negation test to enclose the negations of all individual constraints and then chooses 
the best result to guide the domain splitting during search. A concise description 
of the UCA6 algorithm is presented in Section 4.3. 

For convenience, we define a kind of operator to generalize the idea of a negation 
test, and then give several interesting properties. The generalized operator is called 
the complementary boxing operator, and the corresponding splitting operator is 
called the box splitting operator. 

Definition 3.4 Complementary Boxing Operator, CB. Given a sequence X of n 
real variables associated with domains D. A complementary boxing operator is a 
function CB that takes as input a box x S I" contained in D and a finite set C of 
constraints on X, and that returns a box in I”, denoted by CB(x,C), satisfying the 
following properties: 

(Contractiveness) CB(x,C)Cx, (12) 

(Complementariness) x\CB(x,C)C n (13) 

cgc 



Fig. 4. An example of a complementary boxing (CB) operator applied to a box x and a set {Ci, C 2 } 
of two constraints. 

A box resulting from the application of a complementary boxing operator to 
a bounding box x and a set C of constraints is called a complementary box of C 
within X. The term complementary boxing refers to the process of computing a 
complementary box. The concept of a complementary boxing operator is depicted 
in Figure 4. Additionally, Figure 5 illustrates the outcomes of domain reduction 
operators and complementary boxing operators when applied to the same bounding 
boxes, in some typical situations. 

The complementariness of complementary boxing operators means that the com¬ 
plementary boxing allows isolating certain regions, namely x \ CB(x,C), of which 
the points entirely satisfy all the constraints in C. Especially, if the application of a 
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CB operator = bounding box 



I I = feasible 

DR operator = bounding box DR operator = bounding box 



CB operator 



Fig. 5. Examples of domain reduction (DR) operators and complementary boxing (CB) operators. 

complementary boxing operator to a box and a constraint results in an empty set, 
then the box completely satisfies that constraint. Similarly, if the application of a 
complementary boxing operator to a box with the whole set of constraints results 
in an empty set, then the box is completely feasible. The following theorem states 
this property formally. 

Theorem 3.5. Given a set C of constraints on a sequence of n real variables 
associated with domains T). Suppose x is a box contained in T). If there exists a 
complementary boxing operator CB that maps (x,C) to an empty set (i.e., CB(x, C) = 
%), then C is satisfied with every point in x; that is, x C Hcgc 

A complementary boxing operator can be constructed from a domain reduction 
operator as stated in the following theorem. 

Theorem 3.6. Given a domain reduction operator DR. The function f defined 
as /(x, C) = DR(x, ^C) is a complementary boxing operator. 

Proof. By definition x/ = f{x,C) = DR(x,^C). The contractiveness of domain 
reduction operators implies that xy C x. That is, / enjoys the contractiveness of 
complementary boxing operators. In addition to that, the correctness of domain 
reduction operators implies that 

X n C X/ (14) 

It follows from (14) that, for all a; G x \ xy, we have x ^ x n ^C; thus, x ^ = 

T^\C\c&C ^ ^ HcgC ^ because x G x C P. That is, we have x\x/ C HcgC 

Thus, / enjoys the complementariness of complementary boxing operators. □ 

Theorem 3.7. Given a complementary boxing operator CB. The function f de¬ 
fined as /(x, C) = CB(x, ^C) is a domain reduction operator. 
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Proof. By definition xy = f{x,C) = CB(x, ^C). The contractiveness of comple¬ 
mentary boxing operators implies that xy C x; that is, / enjoys the contractiveness 
of domain reduction operators. In addition to that, the complementariness of com¬ 
plementary boxing operators implies that 

x\x/C^C = I?\ f| C (15) 

cec 

It follows from (15) that, for all a; G x/, we have a; ^ x n thus, x ^'D\ HceC ^ 
and X G HceC ^ because x G Xf C V. That is, we have x\xf C C. This means 
that / enjoys the complementariness of complementary boxing operators. □ 

It follows from Theorem 3.6 and Theorem 3.7 that complementary boxing op¬ 
erators can be constructed from domain reduction operators and vice versa. In 
other words, they are dual to each other. In particular, let C = {Ci,... ,Ck} 
be a set of k constraints. A complementary boxing operator can be constructed 
by CB(x,C) DR(x, ^C) = DR(x, Cq), where Cq is the disjunction of constraints 
-^Ci, ... ,^Ck- In a system that does not accept disjunctive constraints, we can 
relax them by taking the (interval) union of complementary boxes, as stated in 
following theorem. 

Theorem 3.8. Consider a sequence X of n real variables associated with do¬ 
mains V. Let C = {Cl, ... , Ck} be a set of constraints on X and {DRi,..., DR^} 
a set of domain reduction operators for X. Suppose {Cj,..., C{.} is a set of con¬ 
straints on X such that ^Ci = V\ CiCC[ for all i = 1,... ,k. Then the operator 
defined by the following rule is a complementary boxing operator: 

k 

Vx G Io,x C : /(x,C) = nlj DRj(x,C'), (16) 

Proof. The contractiveness of / is obvious because DRi(x, C^,) C x. We now 
prove the complementariness. For every x G x \ DRi(x, C') and i G {1,..., k}: 

X ^ DRi(x, C[) A X n C) 

X ^ X n c' 

X ^ C) 3 ^Ci 
^ xiv\Ci 

^ X G Ci (since x G I?). 

Thus, X G lli^i Therefore, we have 

k k 

x\/(x,C) Cx \ IJdRj(x,C'') C P|0. 

i=l 

This is the complementariness as required. □ 

The negation of a numerical constraint C of the form /(x) o 0 (where o 
is either <, <, >, >, =, or yl) is the constraint /(x) o 0 (where o is either >, 
or =, respectively). In practice, some implementations of domain 
reduction operators only accept constraints that are defined with the relations < 
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and >, but not with the relations < and >. For example, a constraint Ci of the 
form Ci = {f{x) < 0) has the negation of the form -^Ci = {f{x) > 0), which is not 
accepted in some implementations. Fortunately, we can safely use C' = (/(x) > 0) 
in the complementary boxing operator defined by (16) because ^Ci C C' holds. 

The monotonicity of complementary boxing operators is defined similarly to that 
of domain reduction operators (see Definition 3.2). The following theorem gives a 
way to construct a (monotonic) inclusion test. 

Theorem 3.9. Let X be a sequence of n real variables, {DRi,... ,DR„} a set 
of (respectively, monotonic) domain reduction operators and {CBi,..., CB„} a set 
of (respectively, monotonic) complementary boxing operators, where DR^ and CBk 
(k = 1, ... ,n) are defined in k dimensions. Let t be a function that takes as input 
a box X C I" and a finite set C of constraints on a subsequence Y of size k of X, 
and that returns either feasible, infeasible, or unknown such that: 

t(x,C) = infeasible DRfc(x[y],C) = 0; (17) 

t(x.,C) = feasible CBfc(x[y],C) = 0. (18) 

Then t is a (respectively, monotonic) inclusion test. Lf we restrict the codomain 
of T to ezt/ier {infeasible,unknown} or {feasible,unknown}, then the result still 
holds if we use either (17) or (18), respectively, to construct t. 

Proof. It follows directly from Definition 2.4, Definition 2.11, Definition 3.2, 
Theorem 3.3 and Theorem 3.5. □ 

3.3 Domain Splitting Operators 

First, we recall the concept of a bisection, where an interval (i.e., a side) of a box 
is dichotomized into two parts. 

Definition 3.10 Dichotomous Splitting Operator, DS. A dichotomous splitting 
(DS) operator is a function DS : I" ^ 2"° that takes as input a box in I", and 
that returns two (disjoint) boxes in I" resulting from splitting a side of the input 
box into two halves. 

Example 3.11. Consider a box x = (xi,... jXn)"”" C I", where x^ = [x^,Xi[. 
A dichotomous splitting operator DS splits x into two disjoint boxes: x' = 
(xi,... ,Xi_i,x',Xi+i,... ,x„)'’’ and x" = (xi,...,Xi_i, x", x^+i,... ,x„)'’’, where 

= [^i)™cl(x)[ and x" = [mid(x),Xi[. Note that x' n x" = 0; thus, x' n x" = 0. 

Second, we define the concept of a box splitting operator, which splits around a 
complementary box in order to isolate feasible regions w.r.t. a subset of constraints. 

Definition 3.12 Box Splitting Operator, BS. A box splitting (BS) operator is a 
function BS : I" x I" ^ 2'^° which takes as input two boxes such that the former 
contains the latter, and which sequentially splits the outer box along some facets 
of the inner one. The outcome is a sequence of disjoint boxes. 

In fact, the concept of a box splitting operator is a slight generalization of the 
splitting operator proposed in [Van Iwaarden 1996]. The original splitting operator 
gives a way to split a region surrounding a box, provided that this box contains at 
most one optimal solution to a considered optimization problem. 


EPFL Technical Report LIA-REPORT-2006-007, July 2006. 



18 


Xuan-Ha Vu et al. 



CB operator = bounding box 





Fig. 6. Examples of box splitting (BS) and dichotomous splitting (DS) operators. In box splitting, 
all boxes excepted the complementary box are feasible w.r.t. the considered constraints. 

In our algorithm, a box splitting operator that takes as input a domain box 
and a complementary box resulting from the application of a complementary box¬ 
ing operator is applied in combination with a dichotomous splitting operator. The 
dichotomous splitting operator is used when either the complementary boxing oper¬ 
ator produces no reduction or the box splitting operator results in too small boxes. 
Figure 6 illustrates the concept of a box splitting operator for this purpose. 

4. BASIC BRANCH-AND-PRUNE SEARCH ALGORITHMS 

In this section, we first present a generic search in the branch-and-prune framework 
and then present three search algorithms for NCSPs with emphasis on problems 
with continuums of solutions. 

4.1 A Generic Branch-and-Prune Search Algorithm 

We now present a generic search technique, called BnPSearch, for solving NCSPs 
in the branch-and-prune framework, the most common framework for the complete 
solution of NCSPs. The main steps of this generic search technique are described in 
Algorithm 1. Although this technique is not the most general one, it is still capable 
of generalizing the majority of the existing branch-and-prune techniques. 

Taking a CSP as input, the BnPSearch algorithm produces two lists: £v and 
C^. The first list, C\/, is an inner approximation of the solution set. The second 
list, Cs, consists of couples, each consists of a sequence of domains and a set of 
running constraints in these domains. When removing the running constraints, 
will become a boundary approximation of the solution set in association with the 
inner approximation £v- Each couple in Cg constitues a CSP which will need to 
be explored further when reducing the value of the precision e and continuing the 
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solving process. For simplicity, we may omit the running constraints in The 
pruning procedure Prune can be any combination of domain reduction techniques, 
which is not necessarily fixed for all steps of the algorithm. 

Notation 4.1. The notations used in the algorithms in the rest of the paper follow 
the following conventions: 

—The notations B, and B' denote relevant bounding boxes in I"; that is, the 
vectors of domains of the considered NCSPs. 

—The notations C, Ci, C and C denote relevant sets of constraints. 

—The notation CBc denotes a complementary box w.r.t. a constraint, c. 

—The notation and ffl® denote the global lists that accumulate computed boxes 

of inner and boundary union approximations, respectively. 

—The notations DR, CB, and r denote some domain reduction operator, comple¬ 
mentary boxing operator, and inclusion test, respectively. 

4.2 Classic Branch-and-Prune Search Algorithms 

As mentioned in Section 2.3, most complete methods for solving NCSPs integrate 
domain reduction techniques, existence tests, uniqueness tests, exclusion tests, or 
inclusion tests into a bisection search strategy. The most successful solution meth¬ 
ods enhance this process by applying domain reduction techniques such as consis¬ 
tency techniques to the constraint system after each split. This policy is referred 
to as dichotomous maintaining bounds by consistency (DMBC). A generic DMBC 
algorithm is presented in Algorithm 3, where £ is a positive precision (vector) and 
T is a monotonic inclusion test (see Definition 2.11). 


Algorithm 1: BnPSearch - a generic branch-and-prune search algorithm 
Input: a CSP V = (Vo,Oo,Co). 

Output: 'dg. Lists of boxes in inner and boundary union approximations, respectively. 

if PruneCheck(Vo, ©o, Co, e, r, WaitList) then return; -4 On page 19. 

while Waitlist A 0 do 

Get a couple {V, C) from WaitList; t/* Vi g {l,..., fc}: X>i c p, Ci c C. */ 

Split the CSP {Vo,V,C) into sub-CSPs {(Vo, T’i,Ci),..., (Vo, T>fc, Ca,)}; 

for i := 1, . . . ,k do ◄ Do branching, 

if Ci = 0 then 

Cv := Cv U {iDi}; ◄ All points in X>i are solutions 

_ continue for; 

PruneCheck(Vo, Of, Ci, e, r, WaitList); a On page 19. 


Function PruneCheckCV, V, C, e, r, WaitList) 

D := Prune(V, 12, c); 4 Prune the domains in P by using domain reduction techniques, 

if = 0 then return true; 

if all the domains in are not amenable to be split (w.r.t. e) then 
|_ Ce := Ce U {(22, C)}; return true; 
put(WAiTLiST ^ (22,C)); 

return false; 4 The problem has not been solved yet. 
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Algorithm 3: DMBC - an instance of the BnPSearch algorithm 
Input: a bounding box Bq, a constraint set Co, e € Rlfl, a monotonic inclusion test r. 
Output: an inner union approximation a boundary union approximation ffl®. 
if PruneCheckDMBC(Bo, Co, e, T, Waitlist) then return; Page 20 . 

while Waitlist / 0 do 

Get a box B from WaitList; 

(Bi,B 2 ) := Bisect(B); ◄ Vi e {1,2}: Bi c B. 

PruneCheckDMBC(Bi, Co, e, r, WaitList); -4 On page 20 . 

PruneCheckDMBC(B 2 , Co, e, t, WaitList); -4 On page 20 . 


Function PruneCheckDMBCCB, C, e. 

r, WaitList) 

B' :=DR(B,C); 

A Reduce domains. 

if B' = 0 then return true; 

B is infeasible, the problem has been solved. 

if B is canonical or w(B) < e then 



CheckEpsilon(B', C, r); 

< On page 20. 


return true; 


put(WAiTLiST ^ B'); 

Put the current problem into the waiting list. 

return false; 

◄ The problem has not been solved yet. 


Function CheckEpsilonCB, C, r) 

if (Result := r(B,C)) = feasible then 4 Identify the feasibility of B w.r.t. C. 

I ffl^:=ffl^U{B}; B is feasible, store it into the list of inner boxes. 

else if Result = unknown then 


|_ ffl® :=Ea®U{B}; 


B is undiscernible, store it into the list of boundary boxes. 


In most of the known DMBC techniques, the search strategy is to perform split¬ 
ting intervals until the intervals are canonical or their widths are not greater than 
a predetermined precision e. That is, these techniques are able to achieve a pre¬ 
determined precision e (Definition 2.13). In general, the DMBC algorithm cannot 
detect feasible boxes, thus mainly addressing NCSPs with isolated solutions. When 
solving NCSPs with continuums of solutions, we can replace Function PruneCheck- 
DMBC of the DMBC algorithm by Function PruneCheckDMBC+ (on page 21). 
The obtained algorithm is thus called DMBC+. 

The difference (at Line 1) between Function PruneCheckDMBC and Function 
PruneCheckDMBC^ is that the latter resorts to an inclusion test, r' (not t), 
to check if a box is feasible (in other words, it is an inner box). This can also 
be replaced with a complementary boxing operator, CB, checking if it returns an 
empty set. Hence, the DMBC+ algorithm is able to detect feasible boxes, provided 
that the inclusion test r' (or the complementary boxing operator CB) is sufficiently 
efficient. If the inclusion test t' is implemented using an interval form of functions 
as in [Jaulin and Walter 1993; Jaulin et al. 2001], then the DMBC+ will become 
the SIVIA algorithm"* in [Jaulin and Walter 1993; Jaulin et al. 2001]. 

Because of the finite nature of floating-point numbers on computers, it is easy to 

^SIVIA is the abbreviation of set inverter via interval analysis. 
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Function PruneCheckDMBC^ (B, C, e, r, WaitList) 

B^ :=DR(B,C); ◄ Reduce domains, 

if B^ = 0 then return true; B is infeasible, the problem has been solved, 

if B is canonical or w(B) < e then 

CheckEpsilon(B', C, r); -4 On page 20. 

return true; 

1 if r'(B',C) = feasible then -4 This can be optionally replaced with the check CB(B^,C) = 0. 

:= U B is feasible, store it into the list of inner boxes. 

return true; 

put(WAITLlST <— Put the current problem into the waiting list, 

return false; ◄ The problem has not been solved yet. 


prove that both the DMBC and DMBC+ algorithms terminate (i.e., the waiting 
list Waitlist becomes empty) after a finite number of steps and are of the precision 
e w.r.t. the monotonic inclusion test t (see Definition 2.13). 

4.3 New Branch-and-Prune Search Algorithms 

The DMBC and DMBC+ algorithms often generate verbose inner and outer union 
approximations, especially when solving NCSPs with continuums of solutions. The 
first reason is that entirely feasible boxes may be unnecessarily split. The second 
reason is that all constraints are always considered in computations, even when 
some of them are completely satisfied by all points in considered domains. 

Before presenting new search algorithms, we define the following terms. 

Definition 4.2. Given a sequence X ofn real variables {xi, ..., Xn), and a vector 
e = (ei,..., £„)’’" e R” . Consider a box x = (xi,..., Xn)"^ G I” and a finite set C 
of constraints on subsequences of X. A variable Xi is called an active variable in x 
w.r.t. C and e if it occurs in at least one constraint in C, w(xi) = sup(xi) —inf(xi) > 
£i holds and Xi is not canonical. A variable Xi is called an inactive variable in x 
w.r.t. C and e if it is not active in x w.r.t. C and e. 

Inspired by the idea of the ICAb5 algorithm [Benhamou and Goualard 2000], 
we present in Algorithm 7 (on page 22) a simplified version, called UCA5, of the 
UCA6 algorithm [Silaghi et al. 2001], where UCA stands for union-constructing 
approximation. The UCA5 algorithm takes as input an NCSP V = (V,Co,Bo) and 
returns an inner union approximation and a boundary union approximation ffl® 
of the solution set of V. The outer union approximation of the solution set can be 
computed by ffl® := U ffl®. Roughly speaking, the UCA5 algorithm proceeds 
by recursively repeating three main steps/processes: 

(1) Use domain reduction (DR) operators to reduce the current bounding box, which 
plays the role of domains, to a narrower one (see Line 1 in Algorithm 7 and 
Line 5 in Function PruneCheckUCAS). 

(2) Use complementary boxing (CB) operators to search for a complementary box 
w.r.t. some running constraint and the new bounding box obtained at Step 1 
(Line 8 in Function SplitUCAS). During this search, the constraints that make 
empty complementary boxes are removed (Line 10 in Function SplitUCAS). 
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Algorithm 7: UCA5 - a new branch-and-prune search algorithm 


1 

2 

3 


Input: a bounding box Bq, a constraint set Co, e € Rlfl, a monotonic inclusion test r. 
Output: an inner union approximation a boundary union approximation ffl®. 

:= 0; := 0; WAITLIST := 0; ◄ The first two are global lists, 

if PruneCheckUCA5(Bo, Co, e, r, WaitList) then return; -4 On page 22. 

while WaitList / 0 do 

(B,C) := getNext(WAITLlST); -4 Get the next element from the waiting list. 

T = (Splitter, (Bi, ..., Bfc), C, c) := SplitUCA5(B,C); -4 On page 23. 

if T = 0 then continue while; -4 All points in B are solutions, 

for i := 1, . . . ,k do 4 Do branching. 

Ci := C; 

if Splitter = BS und i > 1 then -4 B, does not contain the complementary box. 


Ci := Ci \ {c}; -4 The constraint c is now redundant in Br (Theorem 3.5). 

if Ci = 0 then 4 No running constraints. 

:= U {Bj}; ^ Bi is an inner box. 

_ continue for; 


4 


PruneCheckUCA5(Bi, Ci, e, r, WaitList); 


◄ On page 22. 


Function PruneCheckUCA5(B, C, e, r, WaitList) 

5 B^ := DR(B, C); 4 Prune the domains by using domain reduction operators, 

if B^ = 0 then return true; 4 B is infeasible, the problem has been solved. 

6 if there is no active variable in B' w.r.t. C and e then 4 Definition 4.2. 

CheckEpsilon(B', C, r); 4 On page 20. 

return true; 

7 pLlt(WAITLlST i C))’, 4 Put the current problem into the waiting list, 

return false; 4 The problem has not been solved yet. 


(3) Combine dichotomous splitting (DS) operators with box splitting (BS) operators 
to split the current problem into subproblems (see Line 3 in Algorithm 7 and 
Line 11 to Line 12 in Function SplitUCAS). 

Remark 4.3. In practice, equality constraints usually define surfaces, we then do 
not need to perform the above Step 2 for these constraints. 

The UCA5 algorithm uses a waiting list, WaitList, to store the subproblems 
waiting to be processed further. The elements can be retrieved from, and be put 
to, WaitList by the functions getNext and put. WaitList can be handled as a 
queue or a stack. This allows for the breadth-first search in the former case and 
the depth-first search in the latter case. 

In contrast to the DMBC and DMBC+ algorithms, the UCA5 algorithm restricts 
the DS operators at Line 12 in Function SplitUCAS to dichotomizing a domain of 
a variable only if this variable occurs in at least a running constraint. This avoids 
resulting in a huge number of tiny boxes. The reason is that, in the UCA5 algorithm, 
constraints are removed from consideration whenever empty complementary boxes 
are computed w.r.t. the constraints (see Line 10 in Function SplitUCAS) and that, 
maybe, some variables no longer appear in any running constraints. For simplicity, 
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8 

9 

10 


Function SplitUCA5(B, C) 

foreach C G C do ◄ Search for a constraint c for which complementary boxing results in a reduction. 
CBc i— CB{B, c)i ◄ Enclose the negation of c by using complementary boxing operators. 

if CBc = 0 then 

C C \ {c}; c is redundant in B (Theorem 3.5). 

_ continue for; 

if CBc 7 ^ B then break for; -4 Thus, CBc c B. 

if C = 0 then 4 No running constraints. 

I := U {B}; -4 B is an inner box. 


return 0; 


11 if CBc 7 ^ B then 

L (Bi, ...,B,) ~ BS(B, CBc); Splitter := BS; 4 If bs did not fail, then Bi D CBc. 
if BS failed then Splitter := DS; 

12 if Splitter = DS then (Bi,..., Bfc) := DS(B); 4 Bisect B, fe = 2. 

13 return (Splitter, (Bi, ..., Bfc), C, c); 


the interval (domain) with the greatest width is selected for DS operators. 

For efficiency, the BS operators at Line 11 in Function SplitUCAS split along 
some facet of a complementary box only if this produces sufficiently large boxes; 
the complementary box itself is excepted. This estimation is done by using a 
predetermined parameter, fragmentation ratio. After the splitting phase (Line 13 in 
Function SplitUCAS), if the box splitting operator was chosen and successful (i.e.. 
Splitter = BS), then the first resulting box Bi contains the complementary box 
CBc and the constraint c is always satisfied in all the other boxes because of the 
complementariness of CB operators. 

Function PruneCheckUCAS (on page 22) attempts to apply a DR operator to 
the input subproblem in order to reduce the domains of the subproblem. If this 
returns an empty box, the subproblem has no solutions. Afterwards, the procedure 
at Line 6 in Function PruneCheckUCAS is to check if the input subproblem has 
no active variable. If so, it uses a monotonic inclusion test, called t, to check if the 
subproblem is either infeasible, feasible, or unknown. If r returns infeasible, 
the subproblem is discarded. If t returns unknown, the subproblem is classified as 
undiscernihle w.r.t. e and t. In the other case, every point in the domains of the 
subproblem is a solution. Although the monotonic inclusion test r in our imple¬ 
mentation is a combination of DR and CB operators as described in Theorem 3.9, it 
is however not restricted to this kind of monotonic inclusion test. 

In Algorithm 10, we present a slightly generalized and improved version of the 
UCA6 algorithm - a search technique proposed by Silaghi et al. [2001]. Basically, 
this version is the same as the original version, but it is here improved by changing 
the order of the pruning steps (Function PruneCheckUCAG) such that each box is 
pruned before being put into the waiting list (WaitList). This change reduces the 
number of subproblems in the waiting list because some inconsistent subproblems 
can be discarded sooner than that in the original version in [Silaghi et al. 2001] . This 
version is general enough to be used with the heuristics in [Silaghi et al. 2001] . Those 
heuristics are represented by a generic function, called getSplitType, at Line 11 in 
Function SplitUCAG. Moreover, in this version, we make the stop condition more 
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Algorithm 10: UCA6 - a new branch-and-prune search algorithm 


1 

2 

3 


Input: a bounding box Bq, a constraint set Co, e € Rlfl, a monotonic inclusion test r. 
Output: an inner union approximation a boundary union approximation ffl®. 

:= 0; := 0; WAITLIST := 0; ◄ The first two are global lists, 

if PruneCheckUCA6(Bo, Co, {Bo,... ,Bo}, e, r, WaitList) then return; -4 Page 24. 
while WaitList 0 do ▼/* A set {CBc | c G C} of memorized complementary boxes. */ 

(B,C, (CBc I c G C}) := getNext(WAiTLiST); 

T= (SPLTR,(Bi,...,Bfc),{CBc I cGC}, 6) := SplitUCA6(B,C, {CB,, i c G C}); 

if T = 0 then continue while; -4 All points in B are solutions, 

for i := 1, . . . , fc do 4 Do branching. 

I Ci ;= C; 


if Spltr = BS and i > 1 then 

Ci :— Ci \ {&}; -4 The constraint b is now redundant in the box Br. 

if Ci = 0 then 4 No running constraints. 

:= U {Bi}; 4 Bj is an inner box. 

_ continue for; 


4 


PruneCheckUCA6(Bi, Ci, (CBc | c G Ci}, e, r, WaitList); 4 On page 24. 


Function PruneCheckUCA6(B, C, {CBc | c G C}, e, r, WaitList) 

5 B^ :=DR(B,C); 4 Reduce domains, 

if B^ = 0 then return true; 4 B is infeasible, the problem has been solved. 

6 if there is no active variable in B' w.r.t. C and e then 4 Definition 4.2. 

CheckEpsilon(B', C, r); 4 On page 20. 

return true; 

7 pLlt(WAITLlST <— (B^,C,{CBc I C G C})); 4 Put the current problem into the waiting list, 

return false; 4 The problem has not been solved yet. 


explicit at Line 6 in Function PruneCheckUCAG. It is important (for gaining in 
performance) to emphasize that checking if a domain (i.e., an interval) is not wider 
than the predetermined precision e is only performed on the variables that occur in 
some running constraint, because some constraints have become redundant. This 
detail has been omitted in both [Silaghi et al. 2001] and [Silaghi 2002, Section 5.2.3]. 

In fact, the UCA5 algorithm is a simplification of the UCA6 algorithm. Thus, 
these two algorithms are very much similar and only different at the followings: 

—At Line 1 and Line 4: Function PruneCheckUCAG takes not only the same 
arguments as Function PruneCheckUCAG does but also a list of complementary 
boxes computed at the parent of the current search node. Each box in this list 
is a complementary box w.r.t. a running constraint. At the beginning (Line 1), 
this list consists of boxes that equal to the initial bounding box. 

—At Line 2 and Line 7: Each element in the waiting list WaitList of UCAG 
contains not only a domain box and a set of running constraints but also a list 
of complementary boxes computed at the parent of the current search node. 

—At Line 3: Function SplitUCAG takes not only a box playing the role of domains 
and a set of running constraints but also a list of complementary boxes gotten 
from the waiting list. This function returns not only the used splitting type and 
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Function SplitUCA6(B, C, {CBc | c £ C}) 


8 

9 

10 


foreach c £ C do 

CBc := CB(B n CBc,c); 
if CBc = 0 then C~C\ {c}; 


if C = 0 then 

:= ffl^U{B}; 
return 0; 


11 Splitter := getSplitType(); 
if Splitter = BS then 

CBb := chooseTheBest({CBc | c £ C}); 
(Bi, Bfc) — BS(B, CBfc); 
if BS failed then Splitter := DS; 


c is redundant in B (Theorem 3.5). 

M No running constraints. 
◄ B is an inner box. 

M Get a splitting mode, heuristics can be used. 

M The splitting mode is box splitting. 

M If box splitting did not fail, then Bi D CB^,. 


12 if Splitter = DS then (Bi,..., B*) := DS(B); a Bisect B, fe = 2. 

13 return (Splitter, (Bi, ..., Bj,), {CBc | c £ C}, 6); 


the list of boxes as Function SplitUCAS does but also a new list of complementary 
boxes (because some constraints may have been removed). 

—From Line 8 to Line 10: This is the main difference between the UCA6 algorithm 
and the UCA5 algorithm. While the UCA5 algorithm only finds a complementary 
box that is strictly contained in the bounding box, the UCA6 algorithm computes 
complementary boxes for every constraints and then chooses the best. 

—From Line 11 to Line 13: Since the UCA6 algorithm has just computed com¬ 
plementary boxes for every constraints, it chooses the smallest complementary 
box CBb based on the volume which was computed for the constraint b. The 
constraint b will be used for box splitting operators. In the UCA5 algorithm, the 
first-found complementary box which is strictly contained in the bounding box 
will be used for box splitting operators. Notice that the UCA5 algorithm does 
not need an additional amount of memory to remember the computed comple¬ 
mentary boxes in the waiting list WaitList. 

Of course, both the UCA5 and UCA6 algorithms are instances of the BnPSearch 
algorithm. They compute inner and boundary union approximations, and ffl® 
respectively. These two union approximations are disjoint. Thus, we can obtain an 
outer union approximation by setting ffl® := U ffl®. 

Theorem 4.4. Given a monotonic inclusion test r and a positive precision (vec¬ 
tor) e. Both the UCA5 and UCA6 algorithms terminate and provide inner and 
boundary union approximations, and ffl® respectively, at the precision e with 
respect to the monotonic inclusion test r (see Definition 2.13). 

Proof. No solution is lost because of the correctness of DR operators (Defini¬ 
tion 3.1). Moreover, all points of any boxes in are sound solutions. That is due 
to the complementariness of CB operators (Definition 3.4) and Theorem 3.5. There¬ 
fore, and ffl® are inner and boundary union approximations of the solution set, 
respectively. 

If B is a box in ffl®, then 

—The box B has no active variable w.r.t. e and the running constraints C (at 
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Line 6 of Function PruneCheckUCAS and Function PruneCheckUCA6); 

—Function CheckEpsilon (on page 20) returns unknown (i.e., t(B,C) = unknown). 

If a variable Xi, of which the domain is is in vars(C) (the set of variables of 

C), then the width w(B[a;i]) is not greater than £j. If a variable Xj is not in vars(C), 
we then split the interval B[a;j] of B into intervals whose widths are not greater 
than Sj. Eventually, we obtain a number of e-bounded boxes whose union is B. It 
follows from the properties of a monotonic inclusion test (Definition 2.11) that each 
obtained e-bounded box B' satisfies the property t(B',C) = r(B,C) = unknown, 
since the projections B'[vars(C)] = B[vars(C)]. □ 

Roughly speaking, the larger the union is and the smaller the union ffl® is, 
the more accurate the solution algorithm is. 

5. COMPACTING THE REPRESENTATION OF SOLUTIONS 
5.1 Controlling the Reduction of Small Domains 

When using the UCA5 algorithm or the UCA6 algorithm to solve NCSPs with 
continuums of solutions, we observe that a better alignment of boxes near the 
boundary of the solution set can be obtained by finely controlling the application of 
domain reduction and complementary boxing operators. In particular, if a variable 
Xi of a subproblem is inactive (see Definition 4.2), then the domain reduction or 
complementary boxing should not be enforced on the variable Xj. This is to obtain 
better alignments of contiguous boxes and the computational performance. 

A domain reduction operator (respectively, a complementary boxing operator) 
that only reduce the domains of active variables is called a restricted-dimensional 
domain reduction operator (respectively, a restricted-dimensional complementary 
boxing operator). We denote by DRrd (respectively, CBrd) the restricted-dimensional 
domain reduction operator (respectively, the restricted-dimensional complementary 
boxing operator) obtained from the normal domain reduction operator DR (respec¬ 
tively, the normal complementary boxing operator CB) by enforcing the reduction 
only on active variables. 

When local consistency techniques are enforced in order to obtain the effect of 
domain reduction (i.e., they play the role of domain reduction operators), restricted¬ 
dimensional operators amount to enforcing the local consistency techniques only for 
active variables. The recent algorithms for achieving box consistency, hull consis¬ 
tency and fcB-consistency can be easily modified to adopt the above idea about 
restricted-dimensional operators, for example, by ignoring any procedure involving 
the inactive variables. In case a domain reduction or complementary boxing op¬ 
erator cannot be modified to adopt the idea of reducing only on active variables, 
we can apply it normally and then restore the domains of inactive variables to the 
initial domains. In this case, the gain is only in the (better) alignment of contiguous 
boxes, but not in the performance. 

Fortunately, the implementation of box consistency in a well-known product 
called I LOG Solver [ILOG 2003] supports the above idea about reducing only on 
active variables. It can be done by simply passing only active variables (X) to the 
function IloGenerateBounds when we need to generate narrower bounds on X. 

An illustration of the difference between the effect of a normal domain reduction 
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(a) 


(b) 


Fig. 7 . An example of normal domain reductions and restricted-dimensional domain reductions 
at different levels: (a) all variables {x\ and X2) are considered for the over-reduction; (b) only the 
active variable (3:2) is considered for the reduction. The dark grey regions are inner boxes, the 
light grey regions are undiscernible boxes. 


(dr) operator and the effect of a restricted-dimensional domain reduction (DRrd) 
operator in the solving process is presented in Figure 7. In this example, although 
the normal DR operator produces more accurate output than the DRrd operator does, 
it has to spend much time on making the boundary region narrower than the allowed 
tolerance £i. In practice, this is unnecessary since real world applications mainly 
focus on the inner boxes, the boxes near the boundary are often unsafe for the 
further exploration in the applications, as shown in our experiments (see Section 6). 
Moreover, the number of boxes (15 inner boxes and 8 undiscernible boxes) resulting 
from the application of the DR operator is often higher than the number of boxes 
(3 inner boxes and 8 undiscernible boxes) resulting from the application of the 
DRrd operator. Moreover, the contiguous boxes obtained by using the DRrd operator 
are often aligned; hence, a geometrically compacting technique can work on them 
efficiently to get a concise representation of the solution set, as shown below. 

5.2 Compact Representation of Solutions 

Once the effect of better alignments is obtained, the question is how such a set 
of aligned boxes can be compacted into a smaller set. We propose to use the ex¬ 
treme vertex representation (EVR) of orthogonal polyhedra for this purpose. The 
extreme vertex representation was first proposed by Aguilera and Ayala [1997] for 
the three-dimensional space, and was later extended to the n-dimensional space in 
[Bournez et al. 1999; Bournez and Maler 2000]. The basic idea is that the union of 
disjoint boxes delivered by a box-covering solver defines an orthogonal polyhedron 
for which an improved representation can be generated. An orthogonal polyhedron 
can be naturally represented as the union of disjoint boxes (by enumerating the 
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Fig. 8. Examples of a griddy polyhedron and an orthogonal polyhedron: (a) a griddy polyhedron 
made of the vertex indices of (b) an orthogonal polyhedron. 

boxes and their vertices). That representation is called the disjoint box representa¬ 
tion (DBR) in computational geometry. The EVR is a way of compacting the DBR 
(see [Aguilera 1998] and [Bournez et al. 1999; Bournez and Maler 2000]). Roughly 
speaking, in the EVR, the extreme vertices of an orthogonal polyhedron are iden¬ 
tified and stored in a compact manner such that no information is lost w.r.t. the 
representation of the polyhedra. Moreover, when converting from EVR back to 
DBR, the obtained DBR are often more compact than the initial DBR. 

We recall here basic concepts in the theory of extreme vertex representation. 
The reader can find more details in [Bournez and Maler 2000]. These concepts 
are sufficient to be presented for griddy polyhedra^ because the results on general 
orthogonal polyhedra can be easily obtained from the results on griddy polyhedra by 
mapping the multidimensional array of vertex indices of the orthogonal polyhedra 
to the multidimensional array of vertices of griddy polyhedra (see Figure 8). In 
fact, they do not depend on the orthogonality of the underlying basis. 

For simplicity, polyhedra are assumed to live in X = [0, m[ C (the results 
also hold for X = Rl].). Let G = (0,1,..., m — 1)'^ C be a grid of integer 
points. For every point a; G X, [xj denotes the grid point corresponding to the 
(componentwise) integer part of x. The unit box associated with a grid point 
X = {xi,..., Xd)'^ G 1/ is the box B(a:) = [cci, xi -I- 1[ x • • • x [xd, Xd + 1[. The set 
of all unit boxes is denoted by B. A griddy polyhedron P is the set closure of the 
union of some unit boxes, or can be viewed as a set of grid points. 

Definition 5.1 Color Function. Let P be a griddy polyhedron. The color func¬ 
tion color : X ^ {0,1} is defined as follows: if a: is a grid point then color(a:) = 
1 B(x) C P; otherwise, color(x) = color([a:J). 

We say that a grid point x is black (respectively, white) and that B(a;) is full 
(respectively, empty) when color(x) = 1 (respectively, color(a:) = 0). Figure 9a 
illustrates the color function for griddy polyhedra. Figure 9b illustrates the concept 
of a forward cone based at x G G, which is defined as x^ = {y G G \ x < y}. 

A canonical representation scheme for 2® (or 2^) is a set £ of syntactic objects 

griddy polyhedron is the union of some unit hypercubes with integer-valued vertices (see 
[Bournez et al. 1999]). 
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Fig. 9. A griddy polyhedron: (a) sample colors defined by the color function; (b) the forward cone 

= {y eg \ x <y}. 

such that there is some bijective function 'I' ■. £ —s- 2®; that is, every polyhedron 
has a unique representation. The most simple representation scheme is to explicitly 
enumerate the values of the color function on every grid point; hence, it needs a 
d-dimensional array of bits with m'^ entries. Another simple representation is the 
vertex representation that consists of the set {(a:, color(a;)) | a: is a vertex}. This is 
however still verbose. Hence, it is desired to store only important vertices only. 
The following definition identifies those important vertices. 

Definition 5.2 Extreme Vertex. A grid point x is called an extreme vertex of a 
griddy polyhedron P if the number of black grid points in JV{x) = {xi — l,a;i} x 
• • • X {xd — ^,Xd\ is odd {N{x) is called the neighborhood of x). 

Let 0 denote the exclusive-or (XOR) operation: p (B q = (p A -^q) V (p A -^q). 
The © operation on sets is defined hy A® B = {x \ {x G A) ® (x G B)}. The set of 
extreme vertices (together with the © operation) makes a canonical representation 
of griddy polyhedra as follows [Bournez and Maler 2000, Theorem 1 & 2]. 

Theorem 5.3 Extreme Vertex Representation. For any griddy polyhe¬ 
dron P, there is a unique set V of grid points in Q such that P = 

Moreover, the set V is the set of all extreme vertices of P and this is a canonical 
representation. 

Figure 10 illustrates how the application of the concept of extreme vertex rep¬ 
resentations compacts union approximations. The eight light grey boxes of the 
boundary union approximation in Figure 10a can be concisely and equivalently 
represented by the two light grey boxes in Figure 10b. Theorem 5.3 shows that any 
griddy polyhedron can be canonically represented by the set of its extreme vertices 
(and their colors). 

An effective transformation between DBR and EVR was proposed for low dimen¬ 
sional or small-size (i.e., m is small) polyhedra [Aguilera 1998; Bournez and Maler 
2000]. For example, in the three-dimensional space, the average experimental time 
complexity of converting an EVR to a DBR is far less than quadratic but slightly 
greater than linear in the number of extreme vertices [Aguilera 1998]. Results in 
[Bournez and Maler 2000] also imply that, in a fixed dimension, the time complex¬ 
ity of converting a DBR to an EVR by using the XOR operator is linear in the 
number of boxes in DBR. We then propose to exploit the ideas of these effective 
transformation schemes to produce a compact representation of contiguous aligned 
boxes. This process works as follows: 
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(a) 


(b) 


Fig. 10. The use of extreme vertex representations compacts the boundary union approximation 
in Figure 7b. 

(1) Produce a better alignment of the boxes along the boundaries of constraints. 
This is done by preventing the unnecessary application of reduction operators 
over inactive variables. Figure 7 shows the effect of better alignment obtained 
for a set of nearly aligned boxes of an undiscernible approximation. The original 
set of eight small boxes (Figure 7a) reduces to two groups of four aligned boxes 
(Figure 7b) without altering the predetermined precision. 

(2) The Combination function: The set of aligned boxes in each group, 5i, is 
converted to EVR and then back to DBR to get a set of combined boxes, S 2 
(containing only one box in this case). Due to the properties of EVR, S 2 is 
often more concise than Si. Figure 10 shows how this conversion procedure 
reduces eight boxes to two boxes. 

The above conversion procedure can theoretically be applied in any dimension. 
Due to the efficiency of EVR in low dimension, we however restrict its application 
to very low dimensional small-size regions of the search space in our implementation 
(see Section 5.3). The running time for this conversion is hence near zero. 

5.3 An Improved Branch-and-Prune Search Algorithm 

We present in Algorithm 13 an improved version of the UCA5 and UCA6 algo¬ 
rithms, which is called UCA6+. It also takes as input an NCSP, V = (V,Co,Bo), 
and returns an inner union approximation and a boundary union approxi¬ 
mation ffl® of the solution set of V. Roughly speaking, the UCA6^ algorithm 
uses restricted-dimensional versions of domain reduction and complementary box¬ 
ing operators instead of normal versions in order to produce the effect of better 
alignment (and also to gain in performance), and then uses the conversion between 
the extreme vertex representation and disjoint box representation to get a com¬ 
pact representation of union approximations. The main processes of the UCA6^ 
algorithm are similar to the three main processes of the the UCA5 and UCA6 al- 
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gorithm, but the UCA6^ algorithm uses restricted-dimensional domain reduction 
(DRrd) operators in place of domain reduction (DR) operators (see Line 6 in Func¬ 
tion PruneCheckUCA6+). The UCA6+ algorithm also uses restricted-dimensional 
complementary boxing (CBrd) operators in place of complementary boxing (CB) op¬ 
erators (see Line 13 in Function SplitUCA6+). 

The UCA6^ algorithm does not compute complementary boxes for all running 
constraints as the UCA6 algorithm does. Instead, it allows users to predefine a 
policy to choose a subset C" of C, of which the constraints are enforced with CBrd 
operators (see Line 12 in Function SplitUCA6+). A simple policy is to choose 
either all the constraints of C or a fixed number of constraints in C. A more 


Algorithm 13: UCA6^ - a new branch-and-prune search algorithm 


Input: a box Bq, a constraint set Co, e € R", a monotonic inclusion test r, Dstop- 
Output: an inner union approximation a boundary union approximation ffl®. 

;= 0; := 0; WAITLIST := 0; ◄ The first two are global lists. 

1 if PruneCheckUCA6+(Bo, Cq, 0, e, t, WaitList, Dstop) then return; a On page 31. 
while WaitList / 0 do 

▼ /* A set {CBc I c G C} of complementary boxes that were optionally memorized. */ 

2 (B,C, {CBc I c € C'}) := getNext(WAiTLiST); 4 C' c c, CBc c B. 

3 (SPLTR,(Bi,...,Bfc),{CBc I cGC'},C,b) := SplitUCA6+(B,C,{CBc | c £ C'}); 

if SplitUCAG^ returns 0 then continue while; 

for i := do ◄ Do branching. 

Ci :=C; C'i -C'-, 

if Spltr = BS and i > 1 then 

Ci ~ Ci \ {6}; C'i := C'i \ {6}; 6 is now redundant in the box Bj (Theorem 3.5). 

if Ci = 0 then a No runn ng constraints. 

:= U {Bj}; ^ Bi is an inner box. 

_ continue for; 


foreach c G C'i do if Bi C CBc then C) := C) \ {c}; 
PruneCheckUCA6+(Bi, Ci, {Bi n CBc | c £ C,'}, e, r, WaitList, Dstop); 


Function PruneCheckUCAG^ (B, C, {CBc | c £ C'}, e, r, WaitList, Dstop) 


6 B^ := DRrd(B,C); a Restricted-di mensional domain reduction, 

if B^ = 0 tll©n r©turn trUG^ B is infeasible, the problem has been solved. 

7 if there is no active variable in B' w.r.t. C and e then a Definition 4.2. 

CheckEpsilon(B', C, r); a On page 20. 

return true; 


8 if there are at most Dstop 
(ffl^(B',C), ffl®(B',C)) 


active variables in B' then T/* Resort to another technique. */ 
DimStopSolver(B', C, e, r, DRrd, CBrd); 


T/* Combination(.) does the conversions DBR —EVR —DBR in a Dstop-dimensional space. */ 
:= U Combination(ffl^(B^, C)); ◄ Store in the global list of feasible boxes. 

:= U Combination(ffl^(B^, c)); ◄ Store in the global list of undiscernible boxes. 

r©turn "true^ M The problem has been solved. 


10 foreach c G C' do if B^ C CBc then C' := C' \ {c}; 

11 put(WAITLlST ^ (B^,C, {B' n CBc I C £ ^^}))! the problem into the waiting list, 

return falser M The problem has not been solved yet. 
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Function SplitUCA6+ (B, C, {CBc | c e C'}) 


12 

13 

14 

15 

16 


Choose an arbitrary subset C" ^ C; ◄ C" is a set of constraints to be used with the CBrd operator, 
foreach c € C' U C" do 
if c € C' n C” then 
I CBc:=CBrd(BnCBc,c); 

else if c € C" then < C' . 

I CBc := CBrd(B,c); C':= C'U {c}; 

else < c £ C'. c ^ C". 

|_ CBe — BnCB,; 

if CBc = 0 then C := C \ {c}; c is now redundant in B (Theorem 3.5). 

_ if CBc = 0 or CBc = B then C := C \ {c}; 


17 if C = 0 then ◄ No running constraints. 

:= U {B}; -4 B is an inner box. 

return 0; 


18 

19 

20 


Splitter := g6tSplitTyp6()5 Get a splitting mode, heuristics can be used, 

if Splitter = BS then ◄ The splitting mode is box splitting. 

CBb := chooseTheBest({CBc | c £ C'}); 

(Bi, Bfc) :=BS(B, CB^); M If box splitting did not fail, then Bi D CB 5 . 

if C' = 0 or BS failed then Splitter := DS; 


21 if Splitter = DS then (Bi,..., B*) := DS(B); a Bisect B, fe = 2. 

22 return (Splitter, (Bi, ..., Bfc), {CBc | c e C'}, C, b); 


complicated and dynamic policy based on the pruning efficiency can be used. The 
set of constraints to be considered in the computation of complementary boxes - 
which is done by using CBrd operators (at Line 14 and Line 15) or by intersecting 
with the memorized complementary boxes (at Line 16) - is thus the union C' UC", 
where C is the set of constraints associated with the memorized complementary 
boxes (see Line 2, 3, 5, 11 and 22 in Algorithm 13). Notice that the set C is only 
a subset of C, in general. 

The UCA6+ algorithm uses the same functions getSplitType and chooseTheBest 
as the UCA6 algorithm does (see Line 18 and Line 20 in Function SplitUCA6+). 
The computed complementary boxes can be memorized for improving the com¬ 
plementary boxing of subproblems. However, the memorization should be made 
optional because it may make the computation slow. Unlike the UCA6 algorithm, 
the UCA6^ algorithm only memorizes complementary boxes that do not contain 
the corresponding bounding box (see Line 4 in Algorithm 13 and Line 10 in Function 
PruneCheckUCA6+). 

Function PruneCheckUCA6+ (on page 31) attempts to apply a DRrd operator 
to the input subproblem in order to reduce the domains of this subproblem. If it 
cannot prove that this subproblem is infeasible, it then checks if the subproblem has 
at most Ustop active variables. If the answer is yes, it resorts to a secondary solution 
technique, called DimStopSolver, to solve the current subproblem, provided that 
DimStopSolver provides an output with good alignments. A good candidate for 
DimStopSolver is a search technique with the uniform cell subdivision^ or the 


®A uniform cell subdivision means splitting the domain box into equal e-bounded boxes, called 
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uniform bisection on all variables [Sam-Haroud and Fallings 1996; Lottaz 2000]. 
Variants of the DMBC+ or UCA6 algorithms that use the restricted-dimensional 
operators can also be candidates. 

Given an NCSP, DimStopSolver constructs inner and boundary union approx¬ 
imations of the solution set (see Line 9 in Function PruneCheckUCA6+). These 
two union approximations are naturally represented in DBR (or a bounding-box 
tree). They are converted to EVR and then back to DBR in order to combine each 
group of contiguous aligned boxes into a bigger equivalent box. This conversion 
procedure is performed by the Combination function. 

Theorem 5.4. Given a monotonic inclusion test r and a positive precision (vec¬ 
tor) e. The UCA6+ algorithm terminates and provides inner and boundary union 
approximations, and ffl® respectively, at the precision e with respect to the mono¬ 
tonic inclusion test t (see Definition 2.13). 

Proof. By an argument similar to the proof of Theorem 4.4, we have the fol¬ 
lowing properties: 

(1) and ffl® are inner and boundary union approximations of the solution set, 

respectively; thus, is an outer union approximation of the solution set. 

(2) If not applying the Combination function in Function PruneCheckUCA6+, 
every box B in ffl® can be split into ^-bounded boxes such that each resulting 
box B' satisfies the property: t(B', C) = unknown, where C is the set of running 
constraints in B. 

(3) Since the Combination function does not alter the union of boxes in ffl®, the 
union of boxes in ffl® equals to the union of the above er-bounded boxes. 

This implies what we have to prove. □ 

Notice that all the presented algorithms (DMBC, DMBC+, UCA5, UCA6 and 
UCA6+ ) are complete if the inclusion test t is e-strong for all sufficiently small 
e > 0, since they are all of the precision e w.r.t. r. 

6. EXPERIMENTS 

Since all the above-presented search algorithms should work similarly and equally 
for NCSPs with isolated solutions, in this paper we will only present experiments 
on NCSPs with continuums of solutions. The set of benchmark problems includes 
14 nonlinear problems: PI, P2, P3, P4, FD, G12, H12, F22, LOl, LEI, S06, 
SOS, TD, WP. Their descriptions are given in Appendix A. They are NCSPs with 
continuums of solutions that have been impartially chosen to show different cases 
corresponding to different properties of constraints and solution sets and that can 
be solved efficiently by at least one of the considered search algorithms. 

For the purpose of evaluation, we have implemented five search algorithms 
(DMBC, DMBC+, UCA5, UCA6 and UCA6+) with different options using the 
same data structures and the same domain reduction operators. Our experiments 
discarded DMBC, which is a point-wise approach, as a reasonable candidate for 
solving NCSPs with continuums of solutions because it usually produces a huge 


cells. When this splitting is used, the solver only need to solve subproblems defined on each cell. 

EPFL Technical Report LIA-REPORT-2006-007, July 2006. 



34 


Xuan-Ha Vu et al. 


Table I. The running time results for the search algorithms. The first seven problems are three- 
dimensional while the last seven problems are two-dimensional. 


Algorithm ► 

Splitting type ► 
Memorization ► 

£ 

DMBC+ 

DS 

No 

UCA6 

DS 

Yes 

UCA6+ 

DS 

No 

UCA5 

BS H-DS 

No 

UCA6 

BS -f-DS 

Yes 

UCA6+ 

BS -I- DS 
No 

Ratio 
DMBC + 
UCA6+ 

PI 

lED 

> 24h 

24.94s 

19.34s 

172.74s 

3.76s 

1.03s 

> 83883 

P2 

lED 

> 24h 

187.48s 

95.70s 

6.47s 

7.80s 

0.79s 

> 109367 

P3 

lED 

37724.72s 

61.82s 

25.40s 

8 .86s 

14.23s 

0.63s 

59880 

P4 

lED 

> 24h 

140.25s 

92.89s 

4.96s 

4.51s 

0.96s 

> 90000 

FD 

lED 

505.77s 

183.62s 

101.29s 

48.10s 

59.13s 

31.91s 

15.8 

G12 

lED 

429.82s 

172.52s 

96.40s 

32.72s 

33.59s 

22.23s 

19.3 

H12 

lED 

2161.17s 

889.45s 

267.36s 

280.81s 

273.64s 

99.81s 

21.7 

F22 

IQQI 

5.14s 

3.81s 

3.98s 

3.25s 

3.50s 

2.70s 

1.9 

LOl 

IQQI 

2073.86s 

1082.33s 

660.74s 

49.30s 

51.08s 

7.03s 

295.0 

LEI 

IQQI 

94.04s 

39.79s 

40.89s 

34.35s 

22.05s 

7.32s 

12.8 

S06 

IQQI 

58.58s 

44.29s 

44.84s 

29.78s 

29.10s 

24.34s 

2.4 

SOS 

IQQI 

175.36s 

89.25s 

41.62s 

10.05s 

9.90s 

5.72s 

30.7 

TD 


9.82s 

5.43s 

6.64s 

3.46s 

3.82s 

1.43s 

6.9 

WP 

IQQI 

296.29s 

85.82s 

47.50s 

26.20s 

24.60s 

17.21s 

17.6 


number of boxes, each is £-bounded, in very long running time. The source codes 
of the above search algorithms can be found in the BCS 2.5.2 {box covering solver) 
module, which is downloadable at the official web site of the COCONUT project, 
http://www.mat.univie.ac.at/coconut-environment/. 


Table II. The numbers of boxes in inner union approximations (on the left) and boundary union 
approximations (on the right). 


Prob. 

▼ 

T 

£ 

DMBC+ 

DS 

Memo = No 

UCA6 

DS 

Memo = Yes 

UCA6+ 

DS 

Memo = No 

UCA5 

BS + DS 

Memo = No 

UCA6 

BS -l-DS 
Memo = Yes 

UCA6+ 

BS + DS 

Memo = No 

■a 

■Q 

o 

o 

o 

o 

A 

V 

00 

o 

o 

o 

o 


30601 

222 

15716 

22Sj 

67824 

22^ 

10854 


1253 


■Q 

o 

o 

o 

o 

00 

A 

>730000 


66223 

22^ 

1221 

£ij2] 

23920 


26643 

B2I 

1091 


■Q 

106784 

528757 


2221 

QEQ 

£2^ 


29812 

■£21 

38502 

B2] 

932 

■21 

■Q 

>150000 

>860000 

2^ 

62405 

22U 

31972 

22S 

13988 

22£] 

13423 

B2I 

866 


■Q 

16437 

92681 


92681 

Ifggg 

122] 


65536 

22!! 

70218 

22!! 

35134 


■Q 

17440 

85062 


85062 

2E2I 


2E2] 

59440 

22£I 

60526 

2£2I 

£22] 

mjj 

■Q 

27280 

144296 


144296 


88212 

22&] 

127436 

2£21 

127124 

2]t2j 

£221 

IjgU 

ggj 

1398 

3458 

■gg 

3458 

■jj^ 

^21 

■£21 

2584 

■£2] 

2664 

K2I 

1600 

Qig 

ggj 

65705 

106348 


106348 


£22] 

22!! 

67619 

22&I 

67659 

■2Iq| 

2073 

021 

ggj 

14298 

19688 

22^ 

19688 

22S 

15331 

2^ 

13795 

^££1 

21918 

■£21 

1496 


ggj 

12345 

27756 


27756 


£221 

22S 

30439 

■£21 

26008 

22£j 

|££2I 

gig 

ggj 

26722 

41208 

2jii2] 

41208 

22^ 

32478 

22^ 

27384 

2|q2| 

26624 

|£2I 

11716 

■jg 


2881 

3936 


3936 

■£21 

2915 

2£2! 

4844 

22ii] 

4403 

B2] 

1091 

^■2 

ggj 

22212 

38956 

22£] 

38956 

22!! 

29924 

2E2I 

36433 

■j2£| 

33622 

■£2] 

12!!! 


In the above algorithms, the domain reduction operators (DR and DRrd) have been 
implemented using the function IloGenerateBounds, which is a kind of domain re¬ 
duction and a variant of box consistency, in a well-known commercial product, I LOG 
Solver 6.0 [ILOG 2003]. The complementary boxing operators (CB and CBrd) have 
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Table III. The ratios of the volume of inner approximations to that of outer approximations. 


Algorithm ► 

Splitting type ► 
Memorization ► 

e 

DMBC+ 

DS 

No 

UCA6 

DS 

Yes 

UCA6+ 

DS 

No 

UCA5 

BS -l-DS 

No 

UCA6 

BS + DS 

Yes 

UCA6+ 

BS + DS 

No 

PI 

0.1 

n/a 

0.980 

0.979 

0.997 

0.997 

0.990 

P2 

0.1 

n/a 

0.972 

0.967 

0.996 

0.996 

0.985 

P3 

0.1 

0.710 

0.710 

0.640 

0.956 

0.956 

0.836 

P4 

0.1 

n/a 

0.949 

0.948 

0.997 

0.997 

0.974 

FD 

0.1 

0.984 

0.984 

0.983 

0.992 

0.992 

0.986 

G12 

0.1 

0.874 

0.874 

0.856 

0.924 

0.922 

0.900 

H12 

0.1 

0.885 

0.885 

0.868 

0.938 

0.937 

0.918 

F22 

0.01 

0.968 

0.968 

0.960 

0.977 

0.978 

0.970 

LOl 

0.01 

0.999 

0.999 

0.999 

^1 

^1 

0.999 

LEI 

0.01 

0.997 

0.997 

0.995 

0.999 

0.999 

0.997 

S06 

0.01 

^1 

f«l 

Ril 

^1 

^1 

^1 

SOS 

0.01 

0.999 

0.999 

0.999 

^1 

^1 

^1 

TD 

0.01 

0.996 

0.996 

0.995 

0.998 

0.999 

0.995 

WP 

0.01 

0.999 

0.999 

0.998 

0.999 

0.999 

0.999 


been implemented as in Theorem 3.6. The monotonic inclusion test r has been con¬ 
structed as in Theorem 3.9. The inclusion test r' in Function PruneCheckDMBC+ 
(Page 21) has been implemented using a complementary operator, and thus returns 
either feasible or unknown. For simplicity, the experiments have been taken with 
fixed settings for the new algorithms: fragmentation ratio = 0.25, Ustop = 1- All 
components of the vector e are assumed to be the same. The secondary search 
technique (DimStopSolver) in the UCA6+ algorithm has been implemented as a 
simple combination of a uniform cell subdivision and a monotonic inclusion test. 

The empirical results are presented in Table I, Table II, and Table III. Table I 
shows the running time results of the algorithms (in seconds and hours). Table II 
shows the numbers of boxes in inner and boundary union approximations delivered 
by the algorithms. Table III shows the ratios of the total volume of inner union 
approximations to that of outer union approximations. The term ‘Memorization’ 
(Memo) indicates the memorization of complementary boxes for the next iteration. 
The terms DS and BS -|- DS in Table I, Table II, and Table III indicate the splitting 
policies used in the corresponding search algorithms: 

DS: always dichotomize the largest domain of the domain box. 

BS + DS: attempt to use a box splitting (BS) first; if failed, then proceed with DS. 

Our experiments show that the new algorithms (UCA5, UCA6 and UCA6+) is 
better than the classic algorithm (DMBC+) in all measures. The best gains of the 
new algorithms over the classic one are obtained in case the arities of constraints 
are less than the arity of the problem (e.g., four problems P1-P4). This shows 
how important the reduction of arity of problems is. In most cases, the UCA5 and 
UCA6 algorithms with the option BS-fDS are quite equal in all measures. However, 
the choice of constraints for splitting in UCA6 is far better than that in UCA5 
in the solution of PI. The UCA6+ algorithm with the option BS + DS is always 
better than the others in the running time and the number of boxes, even if it does 
not need the memorization of complementary boxes (thus, less memory is need). 
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The best gains of UCA6+ over the others are obtained when constraint boundaries 
contain a large percentage of nearly axis-parallel regions (e.g., P2 and P3). 

The UCA6+ algorithm with the option BS + DS is slightly less accurate than the 
UCA5 and UCA6 algorithms in the volume measure. However, this situation get 
better when reducing e. Moreover, this is hardly a matter for real world applica¬ 
tions, because no one could ever use all solutions when a very large percentage of 
sound solutions has been found and all the considered algorithms are of precision 
e w.r.t. T (see Definition 2.13). The DMBC+ algorithm and the UCA6 algorithm 
with the option DS produce similar outputs when all constraints in a problem have 
the same set of variables, as happened for all the problems except four problems, 
P1-P4. We observed that the above gains of the new algorithms over the classic 
algorithm, DMBC+, get better when reducing e, especially for hard problems. 

Notice that the arities of constraints in all the above problems, except P1-P4, 
equal to the arities of the problems. In fact, only the experiments on P1-P4 may 
show the full effectiveness of the new algorithms, including the reduction of the arity 
of problem during the solution (we recall that the time and space complexities of 
the algorithms are exponential in the arity of problem). The experiments on the 
other problems do not show the the same improvements as those on P1-P4. This 
reveals that the effect of the arity reduction during the solution in new algorithm is 
an important improvement. Other experiments also show a similar relation among 
the search algorithms using a variant of hull consistency in [Benhamou et al. 1999]. 

7. CONCLUSION 

In this paper, we presented a uniform view on search strategies of branch-and-prune 
methods for solving NCSPs. In this view, we started with a generic branch-and- 
prune algorithm, BnPSearch, and then derived from it two classic algorithms, 
DMBC and DMBC+, for the point-wise and set-covering approaches, respectively. 
As the main contribution of the paper, we proposed three new branch-and-prune 
search algorithms: UCA5, UCA6 and UCA6+. Presenting the new algorithms as 
instances of BnPSearch and extensions of DMBC+ facilitates the comparison of 
them. In particular, we clearly presented the differences among the algorithms. 

Our experiments show that the UCA6+ algorithm (with the option BS -|- DS) is 
the most adaptive search, in time and compactness, among the search algorithms 
for NCSPs with continuums of solutions, while the UCA5 and UCA6 algorithms 
(with the option BS -f- DS) seem to be able to balance between speed and accuracy 
in most cases. They are all far better than the classic branch-and-bound search 
algorithms, DMBC and DMBC+, especially when the arities of constraints are less 
than the arity of the problem. Moveover, the new algorithms often provide a large 
percentage of sound solutions (in the form of a collection or tree of inner boxes) 
when solving NCSPs with continuums of solutions. 

In case the solution set consists of continuums but is highly disconnected, one may 
wish to cluster its union approximations to get grouping/clustering information on 
them. In this case, we propose to use the clustering techniques [Vu et al. 2004] to 
perform a post-processing on the approximations to generate such information. 

We predict that the UCA6+ algorithm will show more speed up and compactness 
if we use higher values, Dgtop =2,3 and use the search technique in [Sam-Haroud 
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and Faltings 1996; Lottaz 2000] in place of DimStopSolver. A direction for further 
research is to explore different combinations of pruning techniques and other tests 
such as existence, uniqueness, exclusion, and inclusion tests to make new branch- 
and-prune search methods, especially for addressing different classes of problems. 
Comparisons with a broader range of search algorithms are also needed. 
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A. NUMERICAL BENCHMARKS 
A.l Problem TD 



^-/ 

L 


Fig. 11. The geometric design of a truss. 


Consider the geometric design problem of a truss depicted in Figure 11. The goal 
is to find the coordinates of the moveable joint {xi,yx) in [0.01,10] x [0.01,10] of 
the node Ni of the truss such that all the following constraints are satisfied: 


F2 

< 

TA, 

A 

< 

TA, 

Fi 

< 

Cl A, 

Fi 

< 

TA, 

Fa 

< 

CaA, 

Fa 

< 

TA, 

lAl 

< 

TA, 

F3 < 0 

=> 

—F3 < C3A, 

Xl 

< 

L, 

yi 

< 

H, 


where 


E 

T 

A 

r 

P 

H 

L 


210 * 10 ® 
235 * 10® 
0.25 
0.5 
400 
6 
10 


Young’s modulus of steel, unit = kN/m^; 

The yield stress of steel, unit = kN/m^; 

The area of cross section of truss members; 

The radius of gyration of the cross section of truss members; 
The loading capacity; 

The height of truss; 

The length of truss; 


and the auxiliary variables are defined as follows: 

tanfi = {H - yi)/{L - xi), 
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tan<:2 

tan ts 
Ri 
Fi 

F2 
Fs 
Fi 
F 5 

Li 

L 3 = + - yif, 

Li = ^x\ 

Ci = T^EKLilrf, 

C3 = TT^EliLsIrf, 

Ci = TT^EKLilrf. 

In fact, this is a two-dimensional problem: the variables are xi and yi. All the 
other variables can be easily eliminated in a preprocessing phase. The reduced 
constraints are however too complicated (in the number of elementary operations) 
to be read, and is hence not listed here. 

A.2 Problem FD 

Consider the design problem of the beam of a railway bridge under cyclic stress. 
The goal is to find {L,qf,Z) G [10,30] x [70,90] x [0.1,10] such that the following 
yield stress and fatigue stress are satisfied:^ 

o- < fy, 

(Te < resistance, 

where 

(Jc = 115000 Yield stress of steel, unit = kN/m^; 

7 = 1.1 The safety factor; 
fy = 460000 Unit = kN/m^; 
years = 200 The number of years to fatigue failure; 

and the auxiliary variables are defined as follows: 

'1.3 ifT<4, 

1.3-0.1(L-4) if4<L<7.5, 

a = 0.95 - 0.008(T - 7.5) if 7.5 < L < 20, 

0.85 - (L - 20)/300 if 20 < T < 50, 

^ 0.75 if L > 50, 

(j) = 0.82-Pl.44/(v^-0.2), 

^The variable Z is scaled up 100 times in unit in comparison to the original version. 


= yi/^1, 

= {H-yi)lxi, 

= PL/H, 

= P/sinti, 

= P/tanti, 

= {Ri - E2)/costs, 

= R\/ COSt2, 

= i?itant2 5 

= J(L-XiY + (H-yiY, 
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Qr 

a 

O-e 

cycles 

CTj' 

resistance 


<7,LV8/(Z/100), 

acr, 

0.05 years, 

(Tc(min{2.5, cycles/2})“^/^, 

(^rh, 


A.3 Problem WP 

This is a two-dimensional simplification of the design model for a kinematic pair 
consisting of a wheel and a pawl. The constraints determine the regions where the 
pawl can touch the wheel without blocking its motion. 

r 20 < + 1/2 < 50; 

I 12y/V(x -12)2 + y2 < 10; 

[ X G [-50,50], y G [0,50]. 

A.4 Problem PI 

Three dimensions; the arities of constraints are less than the arity of problem: 

( 2x^ <3y-{y + 1)°'^ -I- 5; 

I ln(2/3/2 + 2y + l) + 5< z + {z + l/2)0-i; 

I (x -b 1)^-® > 2^/(3 -I- Vz^ + 1); 

[ xG [0,50], yG [0,100], zG [0,50]. 

A.5 Problem P2 

Three dimensions; the arities of constraints are less than the arity of problem: 

{ x^ < y; 
lnj/-b 1 > z; 
xz < 1; 

X G [0,15], y G [1, 200], z G [-10,10]. 

A.6 Problem P3 

P2 added with the fourth constraint whose arity equals to the problem’s arity: 

x^ < y; 

Iny -b 1 > z; 

< xz < 1; 

x^/^ + ln(1.5z -b 1) < y -b 1; 

^ xG [0,15], 2 /G [1,200], zG [0,10]. 

A.7 Problem P4 

Three dimensions; the arities of constraints are less than the arity of problem: 

{ x^-^ -b 1.9 < ln( 2 /^ + y + 1.5); 

ln(2/2 + ^ + 1) < ^ + 2; 

Vx2 -b z2 -b 12x -b 5 < 3 -b (2x -b 3)3; 

X G [0,50], y G [0,100], z G [0,50]. 
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A.8 Problem G12 

Three dimensions; the arities of constraints are equal to the arity of problem: 

( xf + 0.5x2 + 2 (x3 — 3) > 0; 

< xf + X 2 + x"^ < 25; 

[ Xl,X2,X3 e [-8,8]. 

A.9 Problem H12 

Three dimensions; the arities of constraints are equal to the arity of problem: 

{ Xi + X 2 + x'^ < 36; 

(xi — 1)^ + {x2 — 2)^ + > 16; 

+ {x2 — 0 . 4 )^ > 2x3; 

xi,X2,X3 e [-10,10]. 

A.IO Problem F22 

Two dimensions; the intersection of a tricuspoid and a circle: 

( {x^ + 2/^ + 24a: + 36)^ < 64(x + 3)^; 

<1 + 1/2 > 8 ; 

[x,ye [-4,4]. 

A.11 Problem LOl 

Two dimensions; a problem with logarithm and power operations: 

(' (x + 0.1 )y^ > 20 + y/x; 

< ln(VF^ + 13) + 50 > (x + 0.5)1'2; 

[ X G [0,50], y G [0,200]. 

A.12 Problem LEI 

Two dimensions; a problem with logarithm, square root and exponent operations: 

( e^+i/e^^ < 100^/xy + 7 + 30; 

< (x2 — 3x + l)\/y + 2 > Xln(10i/ + 3) + 50; 

[ x,?/ G [0,50]. 

A.13 Problem S06 

Two dimensions; a single constraint whose solution set consists of disconnected 
subsets: 

r 122 //V(x- 12 ) 2 +i/ 2 <io; 

\ X G [-50,50], y G [0,50]. 

A.14 Problem SOS 

Two dimensions; the difference between two circles with interior: 

r 20 < a/x2 +1/2 < 50; 

\ X G [-50,50], y G [0,50]. 


EPFL Technical Report LIA-REPORT-2006-007, July 2006. 



